矩陣乘法優化遞歸式

序:

在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)this


優化遞歸式的特徵

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


本文討論的範疇

  • 形如 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)ip

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


理論基礎

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

已知 :
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)]string

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

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**

相關文章
相關標籤/搜索