矩陣乘法優化遞歸式

序：

在OI比賽中，不少狀況下咱們能夠能經過打表（找規律）或者某些方式發現一個遞歸式。
例如：f(n) = f(n - 1)+f(n - 2)，（斐波那契數列）。

一般狀況下，咱們計算f(n)的時間複雜度就是O(n)（分別計算f(1), f(2) ... f(n - 1)）.
可是當n很大又或者還有其餘處理的複雜度一疊加便會超時。

若是不學習矩陣乘法優化的話，咱們恐怕永遠不會想到計算遞推式還能夠進行優化。
實際上利用矩陣乘法，咱們能夠將O(n)的計算遞歸式的複雜度降至O(logn)。

---

優化遞歸式的特徵

形如f(n) = a1 * f(n - 1) + a2 * f(n - 2) + ... + ak * f(n - k)+c (c爲常數）

---

本文討論的範疇

• 形如 f(n) = f(n-1) + f(n-2) + .. + f(n-k) 的遞歸式(1)

• 形如f(n) = a1 * f(n-1) + a2 * f(n-2) + .. + ak * f(n-k) 的遞歸式(2)

• 形如f(n) = a1 * f(n-1) + a2 * f(n-2) + .. + ak * f(n-k) + c 的遞歸式(3)

實際上理解了最簡單的第一個遞歸式的原理，就很容易理解後面的兩種狀況。每一個前式都是後式的一種特殊狀況。

---

理論基礎

首先給出斐波那契數列求第n項的O(logn)的作法，由它引出原理。

已知 :
f(n) = f(n -1) + f(n - 2); -- (1)
f(n - 1) = f(n - 1); -- (2)
由(1)(2)能夠獲得這樣的一個式子：
[ f(n) ] = [1, 1] * [f(n - 1)]
[f(n - 1)] [1, 0] * [f(n - 2)]

這就核及到矩陣乘法的運算規則。矩陣乘法的計算方式以下所示：

A = X * Y。必須知足row(X) = colom(Y)。
A[i][j] = sum{x[i][k] * x[k][j]}。
row(A) = row(X), colom(A) = colom(Y)..

對於上式，f(n) = f(n - 1) + f(n - 2),f(n - 1) = f(n - 1) + 0，知足斐波那契數列的規則。
咱們設右邊靠左的式子爲A，靠右的爲F。

那麼計算f(n)，咱們只須要計算A^(n - 1) * F。複雜度O(n)。
可是作冪運算，能夠經過快速冪將複雜度從O(n)降到O(logn)，所以總複雜度科技降到O(logn)。

經過這個例子咱們能發現什麼？很顯然的即是A與Y（左式與右二式，一下皆簡稱爲A,X,Y）本質上是一個東西，所以經過迭代，直接計算出第n項。

斐波那契數列是一個最簡單的例子，它也是(1)的典型例子。

廣義斐波那契數列

定義f(n) = a * f(n - 1) + b * f(n - 2).
與以前的區別僅僅在於前面的係數不是1，那麼構造出等式只須要照葫蘆畫瓢便可。

[ f(n) ] = [a, b] [f(n - 1)]
[f(n - 1)] [1, 0] [f(n - 2)],本質沒有變化。
那麼對於f(n) = a1 * f(n - 1) + .... + ak * f(n - k)，咱們只須要構造一個k * k的矩陣。

(*)式即是通常狀況的式子，其實只須要使得第一行知足數列公式，其餘行知足f(i) = f(i)便可。
所以只須要令f[i][i - 1] = 1(下標從0開始)，其餘都是0便可使等式成立。

帶常數與係數的遞歸式

經過以前的經驗咱們不難看出，矩陣A,Y中的元素必定是每一次計算的時候必需的元素。此次多了一個常數c，所以c也要在A,Y中出現。
對此，咱們在須要在最後一行加上c，構形成一個(k + 1) * (k + 1)的矩陣，以下。

觀察一下(**)式，第1和k+1行的最後都爲1，其餘新增的空位全是0，如此便構造完成了。
至於原理，再經過矩陣乘法驗證一下不就行了嗎？

最後給出代碼實現（NOI2012 d1t3）
這道題就是一個裸的(3)式狀況。

/*
About: Matrix Fast Exponentiation to calculate Recursion with Coefficient and Constant
From: NOI2012 d1t3
Description:
    Input: n, p, c, f(0), Mod, mod;
    Recursion: f(n) = (p * f(n-1) + c) % Mod;
    Output: f(n) % mod;
Auther: kongse_qi   
Date:2017/05/19
*/

#include <cstdio>
#include <cstring>

#define read(x) scanf("%d", &x)

typedef long long ll; 

long long mod;

ll Mul (ll a, ll b)
{
    ll ans = 0;
    for(; b; b >>= 1, a = (a + a) % mod) 
    {
        if(b & 1) ans = (ans + a) % mod;
    }
    return ans;
}

struct Matrix{
    ll n, m, x[5][5];
    
    Matrix(){}
    
    Matrix(int a, int b):
        n(a), m(b){}
    
    Matrix operator * (const Matrix& a) 
    {
        Matrix y(n, a.m);
        memset(y.x, 0, sizeof y.x);
        for (int i = 0; i < n; i++)
        {
            for (int k = 0; k < m; k++)
            {
                if (!x[i][k]) continue;
                for (int j = 0; j < a.m; j++)
                {
                    y.x[i][j] = (Mul(x[i][k], a.x[k][j]) + y.x[i][j]) % mod;
                }
            }
        }
        return y;
    }
    
    Matrix Poww(ll times)
    {
        Matrix y = *this, z = y;
        for(--times; times > 0; times >>= 1, z = z * z)
        {
            if(times & 1) y = y * z;
        }
        return y;
    }
};


void Create(Matrix& x, Matrix& f, ll a1, ll p, ll c)
{
    f.n = x.n = x.m = 3;
    f.m = 1;
    memset(x.x, 0, sizeof x.x);
    x.x[0][0] = p;
    x.x[0][2] = x.x[1][0] = x.x[2][2] = 1;
    f.x[0][0] = (Mul(a1, p) + c) % mod; 
    f.x[1][0] = a1;
    f.x[2][0] = c;
    return ;
}

int main()
{
    Matrix x, f;
    ll a1, n, Mod, c, p;
    scanf("%lld%lld%lld%lld%lld%lld", &mod, &p, &c, &a1, &n, &Mod);
    Create(x, f, a1, p, c);
    Matrix y = x.Poww(n - 1) * f; 
    printf("%lld\n", y.x[0][0] % Mod);
    return 0;
}

---

至於一些更簡單的狀況，如下有幾道練習。
luogu P1962 斐波那契數列
luogu P1349 廣義斐波那契數列
Vijos 1067 Warcraft III 守望者的煩惱
NOI2012 隨機數生成器

最後一題值得注意的是，最後幾個點乘法過程會炸出longlong，所以須要使用類快速冪的方法快速求和（積），邊加邊取模。

至此結束。
箜瑟_qi ** 7:58 2017.05.25**