從斐波那契到矩陣快速冪

說起斐波那契數列大家應該都很熟悉,一個簡單的遞推公式

大家應該很容易想出形如這樣的代碼。

int fib(int x) {
    if (x == 1 || x == 2)return 1;
    return fib(x - 2) + fib(x - 1);
}

一個經典的遞歸方法。

但這個代碼的時間複雜度很差,計算到x=40的情況就有點勉強了,因為他其中有太多次重複的計算了。

比如我們輸入x=10,需要計算f(8)與f(9),計算f(9),需要f(8)和f(7),這裡f(8)就重複的計算了兩次,在這重複的遞歸下去重複的計算會越來越多,導致速度緩慢,甚至因為遞歸次數過多導致棧內存溢出(進入函數後會把上一個函數的一系列信息先暫存在棧中,如果遞歸次數過多,導致棧的內存不夠用,就會導致錯誤)。

既然是因為重複計算導致高的時間複雜,我們可以記錄已經計算過的那部分,下次碰到就直接使用,這就是記憶化搜索。

對上面的函數稍加修改可得

//記錄第fb[x]表示第x個斐波那契數列
int fb[100] = {0,1,1};
//為fb數組前3個值賦值
//同fb[0]=0,fb[1]=1,fb[2]=1;

int fib(int x) {
    if (fb[x])return fb[x];//fb[x]非0,說明該值已經計算過,直接返回
    //記憶計算過的值
    fb[x - 1] = fib(x - 1);
    fb[x - 2] = fib(x - 2);
    return fb[x - 1] + fb[x - 2];
}

稍微測試以上的代碼,我們發現速度非常的快了,但發現輸入稍大點的數就會輸出負數,這是因為斐波那契數列往後越來越大,導致超出int,甚至long long的表實範圍。C++沒有java或者python 那樣自帶大數計算,如果需要表現斐波那契數列的後面部分的數值就需要手動寫大數計算,由於斐波那契數列的計算性質我們只要實現大數加法,這倒是不是很麻煩。

大數計算簡單的講就是,使用數組來模擬數字,數組中每一個元素表示數字的一個位。那樣大數加法就是,將每一位相加並進位,就如我們手算加分那樣,乘法同樣是我們從小學的那種乘法豎式計算。

但一般來說不會要我們輸出斐波那契數列的大數,一般來說是對斐波那契數列取模。

稍加修改我們就得出取模的函數

const int mod = 1e9;
int fb[100005] = { 0,1,1 };
int fib(int x) {
    if (fb[x])return fb[x];
    fb[x - 1] = fib(x - 1);
    fb[x - 2] = fib(x - 2);
    return (fb[x - 1] + fb[x - 2])%mod;
}

但我們發現輸入稍大的數就會彈出錯誤,這就是上文中提到的棧溢出,次數過多的遞歸會導致棧溢出,解決方法是,手動模擬遞歸過程,將數據儲存在堆中,或者,換種方式,不使用遞歸計算。模擬遞歸有興趣的可以自行了解,在這就不多展開了。

我們要換種方法,分析遞歸的方法我們發現,遞歸是一種自上而下再從下返回計算的算法。上文的遞歸在自上而下中是沒有進行計算的,只有在自下而上的返回時才有具體的值計算,並且在遞歸中從1~n-1,的斐波那契數列都計算過一遍,那我們能不能直接自下而上的進行計算,顯然可以,我們直接從頭開始循環,一步一步的遞推到結果。

const int mod = 1e9;
int fb[100005] = { 0,1,1 };
int k = 2;
int fib(int x) {
    while (x > k) fb[++k] = (fb[k - 2] + fb[k - 1])%mod;
    return fb[x];
}

觀察遞推公式我們發現,如果沒有多次查詢,那我們一次只要儲存2個數據就可以推出下一個斐波那契數,所以我們還有以下這種寫法

const int mod = 1e9;
int fib(int x) {
    int f[] = { 0,1 };
    int tmp;
    while (--x) {
        tmp = f[1];
        f[1] = (f[1]+f[0])%mod;
        f[0] = tmp;
    }
    return f[1];
}

以上的這幾種算法的時間複雜度已經是O(n)級別了,屬於非常低的一種程度了。

但其實可以進一步優化,使時間複雜度變為O(logn),.斐波那契數列遞推出下一個狀態只需要前兩個數就能推出,也就是說我們有f[]={a,b}這兩個相鄰的斐波那契數,我們需要的是將{a,b}變為{b,a+b}然後重複該操作,也就是上一個代碼的操作。這時就用到一個比較神奇的操作了,我們把數組f[]看成一個2*1的矩陣[a,b],[a,b]到[b,a+b]的變化我們可以用一個矩陣乘法來表示,即[a,b]*A=[b,a+b],顯然當A=時,等式成立。

即[a,b]*=[a*0+b,a+b]= [b,a+b]。

這時我們用數學的方式把一種具體的變化抽象成數學符號。

我們不妨設F(n)=[f(n),f(n+1)]稱它為狀態矩陣(f(n)為第n個斐波那契數列),稱A為轉移矩陣。這時有

F(n+1) = F(n)*A;

即[f(n+1),f(n+2)] = [f(n+1),f(n)+f(n+1)] =[f(n),f(n+1)]*A

根據矩陣乘法的特性易知

F(n)=F(1)*A^(n-1);

這時有的人可能會疑惑了,經行了這麼多的步驟大費周章的把問題數學化,但最後還是要一步一步的去進行矩陣乘法,複雜度甚至比之前的還差。

但是先不要着急,我們這麼大費周章其實是為了這個A,這個不包含可變量的A的連乘,我們有一些方法可以對這個連乘進行優化,也就是對A的次方進行優化。

我們先摒棄A的矩陣屬性,我們把A先化簡成一個普通的數。

這時我們得到一個新的問題,x^y 該如何計算。

大家應該很輕鬆的可以寫出形如這樣的代碼:

int pow(int x, int y) {
    int ans = 1;
    while (y--)ans *= x;
    return ans;
}

代碼複雜度是O(y),我們的目標是優化它,上面代碼的流程是進行累乘,但如果y很大時,一步步累乘繁瑣而緩慢,如果多個多個乘應該會快很多,比如a^7,我們直接寫成a^4+a^2+a 就只要遞推3次就可以得到答案,a^4 可以由a^2自乘得到,a^2也可以由a自乘得到,自乘在代碼中可以使a指數性上升,可以說非常快了,如果運用自乘我們可以將y分割成2^n 的和,如6 分割成 2^2+2^1+2^0。這就與二進制表示不謀而合。我們也可以很輕鬆的用c++/c 中的二進制操作來實現代碼。

int pow(int x, int y) {
    int ans = 1;
    while (y) {
        if (y & 1)ans *= x;
        x *= x;
        y >>= 1;
    }
    return ans;
}

代碼的流程,例如:

輸入a 和23. 23的二進制碼是 10111

先判斷y是否是0,非0進入循環。                                   

y的最低為判斷是否是1,是1就ans乘上x ,10111 最低位是1,所以ans*=x,此時ans=a

然後 a自乘 x=a^2

y右移一位,y=1011

同樣先判斷y是否是0,非0進入循環,取y的最低為判斷是否是1,是1就ans乘上x 1011最低位是1,所以ans*=x,此時ans=a+a^2

直到y=0 跳出循環這時 ans=a+a^2+a^4+a^16=a^23

這就是快速冪,這時複雜度已經達到O(logn)已經是非常快的速度了。按上面的這種方法,類似的我們還可以推出10進制的快速冪,這裡就不做展開了。

 

讓我們回到上面的矩陣冪運算,運用同樣的方法對矩陣使用快速冪,我們就可以得到O(m^3 logn)複雜度 (m是轉移矩陣的大小,m^3 也就是進行一次矩陣乘法的複雜度,)的斐波那契數列算法了。

終於,在講這麼多的前置知識後我們回到我們的正題,矩陣快速冪算法

與上面的快速冪幾乎一樣,只不過乘法換成了矩陣乘法

typedef long long ll;
const ll mod = 1e9 + 7;
const int N = 2;//矩陣大小
int c[N][N];
//狀態轉移
void mul(ll f[N], ll a[N][N]) {
    ll c[N];
    memset(c, 0, sizeof(c));
    for (int j = 0; j < N; j++)
        for (int k = 0; k < N; k++)
            c[j] = (c[j] + f[k] * a[k][j]) % mod;
    memcpy(f, c, sizeof(c));
}
//自乘
void mulself(ll a[N][N]) {
    ll c[N][N];
    memset(c, 0, sizeof(c));
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            for (int k = 0; k < N; k++)
                c[i][j] = (c[i][j] + (a[i][k] * a[k][j])) % mod;
    memcpy(a, c, sizeof(c));
}

//f是狀態矩陣,a是轉移矩陣,n是指數
ll* pow(ll f[N], ll a[N][N],ll n) {

    while (n) {
        if (n & 1)mul(f, a);
        mulself(a);
        n >>= 1;
    }
    return f;
}

例題:

基礎

//vjudge.net/problem/HRBUST-1115   斐波那契數(取模)

//vjudge.net/problem/HRBUST-2170   斐波那契數(大數)

//vjudge.net/problem/POJ-3070  Fibonacci(矩陣快速冪)

提高:

//vjudge.net/problem/HDU-5950     Recursive sequence (acm亞洲賽真題 ,矩陣快速冪運用)

//www.cnblogs.com/komet/p/13770633.html(題解)