矩阵

利用矩阵来求解斐波那契数列的有关问题是ACM题中一个比较常见的题型。例：NYOJ 148(斐波那契数列2)。

有关斐波那契树列的规律详见这里。

(1)、对于n>1，都有f(n)与f(n-1)互质。

(2)、f(n)=f(i)*f(n-i-1)+f(i+1)*f(n-i)。

现在说说怎么利用矩阵来求解斐波那契数列。

我们可以先保存b=f(1),a=f(0)，然后每次设： b'=a+b a'=b。后利用a'和b'一直循环即可。同时我们可以将b a看做一个向量[b a]，前面的操作就可以乘以矩阵：

|1 1|*[b a]=[a+b b]。

|1 0|

也就是说，如果我们要求第100个fibonacci数，只需要将矩阵[1 0]乘上 

1 1 

1 0 

的一百次方，再取出第二项即可。就像题目中的描述一样：设F₀ = 0，F₁ = 1，那么有：

.

这个矩阵的n次方的形式是：

F(n+1)  F(n)

F(n)   F(n-1)，

现在的问题是如何求出这个矩阵的n次方是多少？

可以使用的方法有两种：

(1)、可以利用二分的思想：假设要求矩阵的N次方为A(N)，设i=N/2若N%2==1，则 A(N)=A(i)*A(i)*A(1)若N%2==0，则A(N)=A(i)*A(i)。

(2)、利用二进制的思想：将N拆为二进制数，譬如13=(1101)2那么 A^13= A^8 * A^4 * A^1。这个在求快速幂模的时候经常用。具体见代码，可做模板，写的比较的乱：

[cpp] view plaincopyprint?

1. #include<iostream>  
2. #include<cstring>  
3. using namespace std;  
4. const int MAX=10;  
5. #define __int64 long long  
6. #define Bit(n) 1<<n  
7. #define CLR(arr,val) memset(arr,val,sizeof(arr))  
8. class Matrix{  
9. public:  
10.     Matrix(int r,int c):row(r),col(c){}  
11.     void Init()    
12.     {   CLR(map,0);  
13.         map[0][0]=map[0][1]=map[1][0]=1;  
14.     }  
15.     void Unit() //初始化为单位矩阵   
16.     {   CLR(map,0);  
17.         for(int i=0;i<row;i++)  
18.             map[i][i]=1;  
19.     }  
20.     int Result() const{return map[0][1]%10000;}  
21.     friend Matrix operator*(const Matrix& ,const Matrix&);  
22.     int Pow(int);  
23. private:  
24.     __int64 map[MAX][MAX];      
25.     int row,col;         
26. };  
27. Matrix operator*(const Matrix& M1,const Matrix& M2) //矩阵相乘模板  
28. {   Matrix M(M1.row,M2.col); //相乘之后矩阵的行和列会变化     
29.     for(int i=0;i<M1.row;i++)  
30.         for(int j=0;j<M2.col;j++)  
31.         {   M.map[i][j]=0;  
32.             for(int k=0;k<M1.col;k++)  
33.                 M.map[i][j]+=M1.map[i][k]*M2.map[k][j];   
34.             M.map[i][j]%=10000;  
35.         }    
36.     return M;      
37. }  
38. Matrix M(2,2);  
39. int Matrix::Pow(int n) //矩阵快速幂   
40. {   Matrix temp(2,2);   
41.     temp.Init();   
42.     for(int i=0;Bit(i)<=n;i++) //利用二进制的思想求解   
43.     {   if(Bit(i)&n) M=M*temp;   
44.         temp=temp*temp;       
45.     }  
46.     return M.Result();      
47. }  
48. int main()  
49. {   __int64 num;  
50.     while(cin>>num,num!=-1)  
51.     {   M.Unit();  
52.         cout<<M.Pow(num)<<endl;  
53.     }  
54.     return 0;  
55. }  

在这里顺便贴下斐波那契数列的通项公式：

 我们可以利用这个公式来简化很多运算：例HDU 2855，意思是输入n,m，求%m的结果。

具体推导过程如下：

 

 

 

 

 

 

 

 

 

 

 

 上个题目中求F(n)%m的值，这个是求F(2n)%m的值，你懂的~代码就不多写了~

调试了几次10^9这个数就是过不了，矩阵快速幂的模板改下就过了~具体改成下面的：

[cpp] view plaincopyprint?

1. int Matrix::Pow(int n) //矩阵快速幂   
2. {   Matrix temp(2,2);   
3.     temp.Init();   
4.     while(n) //利用二进制的思想求解   
5.     {   if(n&1) M=M*temp;   
6.         temp=temp*temp;       
7.         n>>=1;  
8.     }  
9.     return M.Result();      
10. }  

例题2：NYOJ 507(斐波那契数列)，只要写出那个递推关系式即可~其余的直接套模板：

S(n)=f(0)²+f(1)²+f(2)²+...+f(n)²；

S(n-1)=f(0)²+f(1)²+f(2)²+...+f(n-1)²。 

所以：S(n)-S(n-1)=f(n)²=x²f(n-1)²+y²f(n-2)²+2xyf(n-1)f(n-2)。

所以可以构造下面的矩阵：

这个时候就简单了，具体见代码：

[cpp] view plaincopyprint?

1. #include<iostream>  
2. #include<cstring>  
3. #include<cstdio>  
4. using namespace std;  
5. const int MAX=10;  
6. #define __int64 long long   
7. #define CLR(arr,val) memset(arr,val,sizeof(arr))  
8. int x,y,num;  
9. class Matrix{  
10. public:  
11.     Matrix(int r,int c):row(r),col(c){}  
12.     void Init()    
13.     {   CLR(map,0);  
14.         x%=10007,y%=10007;//一定要把这个条件加上~~  
15.         map[0][0]=map[0][1]=map[2][1]=1;   
16.         map[1][1]=x*x%10007;  
17.         map[1][2]=y*y%10007;  
18.         map[1][3]=(2*x*y)%10007;  
19.         map[3][1]=x%10007,map[3][3]=y%10007;   
20.     }  
21.     void Unit() //初始化为单位矩阵   
22.     {   CLR(map,0);  
23.         for(int i=0;i<row;i++)  
24.             map[i][i]=1;  
25.     }  
26.     void Clear()  
27.     {   for(int i=0;i<4;i++)  
28.             map[i][0]=1;  
29.     }  
30.     int Result() const {return map[0][0]%10007;}  
31.     friend Matrix operator*(const Matrix& ,const Matrix&);  
32.     friend ostream& operator<<(ostream& ,const Matrix&);   
33.     int Pow(int);  
34. private:  
35.     __int64 map[MAX][MAX];      
36.     int row,col;         
37. };  
38. Matrix operator*(const Matrix& M1,const Matrix& M2)   
39. {   Matrix M(M1.row,M2.col);  
40.     for(int i=0;i<M1.row;i++)  
41.         for(int j=0;j<M2.col;j++)  
42.         {   M.map[i][j]=0;  
43.             for(int k=0;k<M1.col;k++)  
44.                 M.map[i][j]+=M1.map[i][k]*M2.map[k][j];   
45.             M.map[i][j]%=10007;  
46.         }    
47.     return M;      
48. }  
49. Matrix M(4,4);  
50. int Matrix::Pow(int n)   
51. {   Matrix temp(4,4),Start(4,1);  
52.     temp.Init();  
53.     Start.Clear();  
54.     //cout<<temp;   
55.     while(n)  
56.     {   if(n&1) M=M*temp;  
57.         temp=temp*temp;  
58.         n>>=1;  
59.     }  
60.     Start=M*Start;  
61.     return Start.Result();      
62. }  
63. ostream& operator<<(ostream &cout,const Matrix& M)  
64. {   for(int i=0;i<M.row;i++)  
65.     {   for(int j=0;j<M.col;j++)  
66.             cout<<M.map[i][j]<<" ";  
67.         cout<<endl;      
68.     }  
69.     return cout;  
70. }  
71. int main()  
72. {   while(scanf("%d%d%d",&num,&x,&y)!=EOF)  
73.     {   M.Unit();  
74.         cout<<M.Pow(num)<<endl;  
75.     }  
76.     return 0;  
77. }