從斐波那契到矩陣快速冪

_comet發表於2020-10-09

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

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

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;
}

例題:

基礎

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

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

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

提高:

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

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

 

相關文章