1. 程式人生 > >矩陣快速冪求斐波那契數列(初學整理)

矩陣快速冪求斐波那契數列(初學整理)

  參考文章:

         感謝以上兩位大神讓我明白瞭如何用矩陣快速冪求斐波那契數列。

  題外話:線代真的要好好學,這學期剛好學線代,初期還是蠻認真的,可是後來就勉強能聽懂,然後慢慢~~。萬惡的線代,上學期去湘潭大學參加中南地區程式設計邀請賽第一題就是一個矩陣,三個人都沒有學過線代,可那題的通過率還是蠻高的,最終只做了一題打鐵的佳績!

         感覺學起來還是蠻簡單的,蠻好理解的,以下就我的理解借用大神們一點想法來把矩陣快速冪整理總結!

第一部分:矩陣的基礎知識

1.結合性 (AB)C=A(BC).

2.對加法的分配性 (A+B)C=AC+BCC(A+B=CA+CB 

3.對數乘的結合性 k(AB=kA)B =A(kB).

4.關於轉置 (AB)'=B'A'

一個矩陣就是一個二維陣列,為了方便宣告多個矩陣,我們一般會將矩陣封裝一個類或定義一個矩陣的結構體,我採用的是後者。(弱雞的我也直只會用結構體實現)

第二部分:矩陣相乘

若A為n×k矩陣,B為k×m矩陣,則它們的乘積AB(有時記做A·B)將是一個n×m矩陣。前一個矩陣的列數應該等於後一個矩陣的行數,得出的矩陣行數等於前一個矩陣的行數,列數等於後一個矩陣的行數。

其乘積矩陣AB的第i行第j列的元素為:

舉例:A、B均為3*3的矩陣:C=A*B,下面的程式碼會涉及到兩種運算順序,第一種就是直接一步到位求,第二種就是每次求一列,比如第一次,C00+=a00*b00,C01+=a00*b01……第二次C00+=a00*b10,C01+=a01*b11……以此類推。。。

C00 = a00*b00 + a01*b10 + a02*b20
C01 = a00*b01 + a01*b11 + a02*b21 
C02 = a00*b02 + a01*b12 + a02*b22
C10 = a10*b00 + a11*b10 + a12*b20
C11 = a10*b00 + a11*b11 + a12*b21
C12 = a10*b02 + a11*b12 + a12*b22
C20 = a20*b00 + a21*b10 + a22*b20
C21 = a20*b01 + a21*b11 + a22*b21
C22 = a20*b02 + a21*b12 + a22*b22
C00 = a00*b00 + a01*b10 + a02*b20
C01 = a00*b01 + a01*b11 + a02*b21 

C02 = a00*b02 + a01*b12 + a02*b22

下面先來實現一個矩陣相乘的函式吧。

const int MOD=10000;
struct mat
{
    int a[2][2];//這裡資料範圍就用小的示範
};
mat mat_mul(mat x,mat y)//實現兩個矩陣相乘,返回的還是一個矩陣。
{
    mat res;//用來表示得到的新的矩陣;
    memset(res.a,0,sizeof(res.a));
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        for(int k=0;k<2;k++)
    {
        res.a[i][j]+=x.a[i][k]*y.a[k][j];
        res.a[i][j]%=MOD;//這一步看題目具體需要了
    }
    return res;
}

學了現代的話這個是很好理解的(個人認為)。

第三部分:矩陣快速冪   //其實和普通快速冪類似,只不過這裡需要得到的是一個矩陣

神馬是冪?【很多時候會被高大上的名字嚇到。。。導致學習效率降低。。。其實沒辣麼可怕,很簡單!!!】

冪又稱乘方。表示一個數字乘若干次的形式,如n個a相乘的冪為a^n ,或稱a^n為a的n次冪。a稱為冪的底數,n稱為冪的指數。——引自.度娘百科

這類題,指數都是很大很大很大很大很大很大很大的。。。霸王硬上弓的話,很容易超時的 T_T 。。。所以得快速冪→_→

學過之後發現,其實矩陣快速冪 的核心思想跟 以前學過的快速冪取模非常非常相似,只是矩陣乘法需要另外寫個函式,就是上面那個程式碼。。。

快速冪的思路就是:

設A為矩陣,求A的N次方,N很大,1000000左右吧。。。

先看小一點的,A的9次方

A^9

= A*A*A*A*A*A*A*A*A  【一個一個乘,要乘9次】

= A*(A*A)*(A*A)*(A*A)*(A*A)【保持格式的上下統一,所以加上這句】

 = A*(A^2)^4 【A平方後,再四次方,還要乘上剩下的一個A,要乘6次】

= A*((A^2)^2)^2【A平方後,再平方,再平方,還要乘上剩下的一個A,要乘4次】

也算是一種二分思想的應用吧,1000000次冪,暴力要乘1000000次,快速冪就只要(log2底1000000的對數) 次,大約20次。。。這。。。我沒錯吧。。。

單位矩陣: n*n的矩陣 mat ( i , i )=1; 任何一個矩陣乘以單位矩陣就是它本身 n*單位矩陣=n, 可以把單位矩陣等價為整數1。(單位矩陣用在矩陣快速冪中)

例如下圖就是一個7*7的單位矩陣:



下面來實現一個矩陣快速冪:

int pow(int n)//還是小範圍資料來說吧,要不然返回值的型別自己定義
{
    mat c,res;
    memset(res.a,0,sizeof(res.a));
    c.a[0][0]=1;//給矩陣賦初值
    c.a[0][1]=1;
    c.a[1][0]=1;
    c.a[1][1]=0;
    for(int i=0;i<n;i++) res.a[i][i]=1;//單位矩陣;
    while(n)
    {
        if(n&1) res=mat_mul(res,c);//這裡看就要用到上面的矩陣相乘了;
        c=mat_mul(c,c);
        n=n>>1;
    }
    return res.a[0][1];
}//時間複雜度log(n)

      但是矩陣如何與斐波那契聯絡在一起呢???

找了很多部落格,看了第二位大神的部落格才理解。

對於矩陣乘法與遞推式之間的關係:

如:在斐波那契數列之中

f[i] = 1*f[i-1]+1*f[i-2]  f[i-1] = 1*f[i-1] + 0*f[i-2];

所以

就這兩幅圖完美詮釋了斐波那契數列如何用矩陣來實現。

  給出了矩陣相乘的定義,要你求出斐波那契的第n項對1e4取餘。

  程式碼一:在網上看到的很簡潔的程式碼

#include <iostream>
#include <cstddef>
#include <cstring>
#include <vector>
using namespace std;
typedef long long ll;
const int mod=10000;
typedef vector<ll> vec;
typedef vector<vec> mat;
mat mul(mat &a,mat &b)//表示不會這樣用,,,,
{
    mat c(a.size(),vec(b[0].size()));
    for(int i=0; i<2; i++)
    {
        for(int j=0; j<2; j++)
        {
            for(int k=0; k<2; k++)
            {
                c[i][j]+=a[i][k]*b[k][j];
                c[i][j]%=mod;
            }
        }
    }
    return c;
}
mat pow(mat a,ll n)
{
    mat res(a.size(),vec(a.size()));
    for(int i=0; i<a.size(); i++)
        res[i][i]=1;//單位矩陣;
    while(n)
    {
        if(n&1) res=mul(res,a);
        a=mul(a,a);
        n/=2;
    }
    return res;
}
ll solve(ll n)
{
    mat a(2,vec(2));
    a[0][0]=1;
    a[0][1]=1;
    a[1][0]=1;
    a[1][1]=0;
    a=pow(a,n);
    return a[0][1];//也可以是a[1][0];
}
int main()
{
    ll n;
    while(cin>>n&&n!=-1)
    {
        cout<<solve(n)<<endl;
    }
    return 0;
}

程式碼二:自己寫的程式碼
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD=10000;
struct mat
{
    ll a[2][2];
};
mat mat_mul(mat x,mat y)
{
    mat res;
    memset(res.a,0,sizeof(res.a));
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        for(int k=0;k<2;k++)
        res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j])%MOD;
    return res;
}
void mat_pow(int n)
{
    mat c,res;
    c.a[0][0]=c.a[0][1]=c.a[1][0]=1;
    c.a[1][1]=0;
    memset(res.a,0,sizeof(res.a));
    for(int i=0;i<2;i++) res.a[i][i]=1;
    while(n)
    {
        if(n&1) res=mat_mul(res,c);
        c=mat_mul(c,c);
        n=n>>1;
    }
    printf("%I64d\n",res.a[0][1]);
}
int main()
{
    int n;
    while(~scanf("%d",&n)&&n!=-1)
    {
        mat_pow(n);
    }
    return 0;
}

  感覺不是很難,裸模板的話沒什麼問題,但如果靈活度高一點我就不會構造矩陣了。還是繼續加油吧!

  再次感謝大神們的部落格。