矩陣快速冪專題（一）
最近閒來無事，準備集中精力刷一波數論與圖論。矩陣快速冪是數論裡面的重要組成部分，值得我好好學習一下。因為題目比較多，分析也比較多，所以將此專題分成幾個部分。做完這一專題，可能會暫時轉向圖論部分，然後等我組合數學學得差不多了，再回過頭來繼續做數論題。
矩陣快速冪演算法的核心思想是將問題建模轉化為數學模型（有一些簡單題目是裸的矩陣模型，但是大部分難題就是難在要構造矩陣，用矩陣方法解決問題），推倒遞推式，構造計算矩陣，用快速冪的思想求解矩陣A的n次方取mod，從而得到矩陣裡面你需要的資料。
矩陣快速冪問題有某些特徵，是類似快速冪的形式。首先，從題目中是一定可以得到一個遞推式的，而且這個遞推式不是簡單的An與An-1的關係，而是關係式中帶有An-2或之後的項，還有的情況甚至帶有常數項和多個數列的組合（加減乘數）；其次，資料量比較大，不可能用打表方法做出來；最後，有些題目可能與矩陣一點沒有關係，但是通過規律分析，能夠轉化為已知遞推式求解某一項（取模）的。
這裡給出模版： 我覺得我的矩陣模版不是特別好，但是也不是特別麻煩，所以請讀者自行選擇是否採用
#define maxn 10+5; typedef long long ll; const ll mod = 1e9+7; #define clr(x,y) memset(x,y,sizeof(x)) struct matrix { int n; ll maze[maxn][maxn]; void init(int n) { this->n=n; clr(maze,0); } matrix operator *(const matrix& rhs) { matrix ans; ans.init(n); for(int i=0;i<n;i++) for(int j=0;j<n;j++) for(int k=0;k<n;k++) ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod; return ans; } }; matrix qlow(matrix a,int n) { matrix ans; ans.init(a.n); for(int i=0;i<a.n;i++)ans.maze[i][i]=1; while(n) { if(n&1)ans=ans*a; a=a*a; n>>=1; } return ans; }
第一題 hdu-1757
分析：題目直接給出了遞推式（帶有f（x-1）到f（x-10）），而且資料量非常大需要取模(這裡並不要管a0~a9的值，若不為01那麼只需多考慮一步取模)，這是簡單的矩陣快速冪裸題。
第二題 hdu-1575#include<set> #include<map> #include<cmath> #include<stack> #include<queue> #include<vector> #include<iostream> #include<cstdio> #include<cstring> #include<string> #include<algorithm> #include<cctype> #define maxn 6+5 #define clr(x,y) memset(x,y,sizeof(x)) using namespace std; const int inf = 0x3f3f3f3f; typedef long long ll; const double pi = acos( -1 ); //const ll mod = 1e9+7; int a[10]; ll k,m; ll qmul(ll a,ll b) { ll ans=0; while(b) { if(b&1)ans=(ans+a)%m; a=(a+a)%m; b>>=1; } return ans; } struct matrix { int n; ll maze[maxn][maxn]; void init(int n) { this->n=n; clr(maze,0); } matrix operator * (matrix & rhs) { matrix ans; ans.init(n); for(int i=0;i<n;i++) for(int j=0;j<n;j++) for(int k=0;k<n;k++) ans.maze[i][j]=(ans.maze[i][j]+qmul(maze[i][k],rhs.maze[k][j]))%m; return ans; } }; matrix qlow(matrix a,ll n) { matrix ans; ans.init(a.n); for(int i=0;i<ans.n;i++)ans.maze[i][i]=1; while(n) { if(n&1)ans=ans*a; a=a*a; n>>=1; } return ans; } int main() { //freopen("d:\\acm\\in.in","r",stdin); while(~scanf("%lld %lld",&k,&m)) { for(int i=0;i<10;i++) scanf("%d",&a[i]); if(k<10) { printf("%lld\n",k%m); continue; } matrix ans; ans.init(10); for(int i=0;i<10;i++)ans.maze[0][i]=i; matrix ant; ant.init(10); for(int i=0;i<9;i++)ant.maze[i+1][i]=1; for(int i=0;i<10;i++)ant.maze[i][9]=a[9-i]; matrix tmp; tmp=qlow(ant,k-9); ans=ans*tmp; printf("%lld\n",ans.maze[0][9]); } return 0; }
分析：這道題就是非常直白的矩陣快速冪，將主對角線上的數加起來取模-。-
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 6+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 9973;//1e9+7;
ll qmul(ll a,ll b)
{
ll ans=0;
while(b)
{
if(b&1)ans=(ans+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return ans;
}
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (matrix& rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=(ans*a);
a=(a*a);
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
int t;
scanf("%d",&t);
while(t--)
{
int n,k;
scanf("%d %d",&n,&k);
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
scanf("%d",&ans.maze[i][j]);
ans=qlow(ans,k);
ll anw=0;
for(int i=0;i<n;i++)
anw=(anw+ans.maze[i][i])%mod;
printf("%lld\n",anw);
}
return 0;
}分析：我這裡只提出我的方法，其他網上大神的方法請到別的部落格上查詢。
對於L較大的情況（L大於2，至於為什麼是2，你沒有看見fmf和fff都是三個嗎？），我考慮對於在E-queue裡面的所有合法佇列（-。-，其實和佇列沒有任何關係）。他們都是以mm，fm，mf，ff之一結尾的，我就分別將以他們結尾的並且總長度為n的個數記為mm[n]，fm[n]，mf[n]，ff[n]。
考慮mm[n]後面加m是mm[n+1]的一部分，加f是mf[n+1]的一部分；
fm[n]加f屬於O-queue，捨去，加m是mm[n+1]的一部分；
依次類推。。。。。。
得到遞推式mm[n+1]=mm[n]+fm[n]; mf[n+1]=mm[n+1]; fm[n+1]=ff[n]+mf[n]; ff[n+1]=mf[n+1];
構造矩陣
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;
int m;
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix& rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%m;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
int l;
while(~scanf("%d %d",&l,&m))
{
if(l<=2)
{
if(l==0)puts("0");
else if(l==1)printf("%d\n",2%m);
else printf("%d\n",4%m);
continue;
}
matrix ans;
ans.init(4);
for(int i=0;i<4;i++)ans.maze[0][i]=1;
matrix ant;
ant.init(4);
ant.maze[0][1]=ant.maze[1][3]=ant.maze[2][0]=ant.maze[2][1]=ant.maze[3][2]=ant.maze[3][3]=1;
ant=qlow(ant,l-2);
ans=ans*ant;
int anw=0;
for(int i=0;i<4;i++)anw=(anw+ans.maze[0][i])%m;
printf("%d\n",anw);
}
return 0;
}第四題 hdu-2256
分析：一開始會以為與矩陣無關，但是看見根號還要取模就會發現滿滿都是套路。-。-雖然根號2與根號3是兩個根號，但是還有個2次方，先將2次方乘進去，得到5+2sqrt(6)，發現根號減少到一個，而且對於5+2sqrt(6)的n次方，一定是一個常數項和一個帶根號6的項，而且互相之間有遞推關係，這就把問題引向矩陣上了。
但是這裡不能直接取模，因為Bn還有一部分是sqrt(6)，直接取模顯然是錯的
這樣構造的矩陣和最後要求的資料都有了 -。-
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1024;//1e9+7;
struct matrix
{
int n;
ll maze[2][2];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix& rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
int n;
scanf("%d",&n);
matrix ans;
ans.init(2);
ans.maze[0][0]=1;
matrix ant;
ant.init(2);
ant.maze[0][0]=ant.maze[1][1]=5;
ant.maze[0][1]=2;
ant.maze[1][0]=12;
ant=qlow(ant,n);
ans=ans*ant;
printf("%d\n",(2*ans.maze[0][0]-1)%mod);
}
return 0;
}第五題 codeforces-185A
分析：做完前面的題，這道題反而簡單了。不妨考慮當n年後（操作n次），正三角的個數是An，倒三角的個數是Bn。那麼繼續再操作一次，An個正三角變成了3An個正三角和An個倒三角，Bn個倒三角變成Bn個正三角和3Bn個倒三角。所以有遞推關係 An+1=3An+Bn Bn+1=An+3Bn
構造矩陣
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;
ll mul(ll a,ll b)
{
ll ans=0;
while(b)
{
if(b&1)ans=(ans+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return ans;
}
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix& rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+mul(maze[i][k],rhs.maze[k][j]))%mod;
return ans;
}
};
matrix qlow(matrix a,ll n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
ll n;
scanf("%I64d",&n);
matrix ans;
ans.init(2);
ans.maze[0][0]=1;
matrix ant;
ant.init(2);
ant.maze[0][0]=ant.maze[1][1]=3;
ant.maze[1][0]=ant.maze[0][1]=1;
ant=qlow(ant,n);
ans=ans*ant;
printf("%I64d\n",ans.maze[0][0]);
return 0;
}分析：題目講的好像是九連環（好像又不是，英語不好），-。- 可憐我並沒有玩過
這道題非常像漢諾塔，但是題目只是求最小步數，而不需要求出具體的步驟。我假設對於長度為i的中國環，我最少需要f(i)步（1<=i<=n）。那麼對於長度為n+1的中國環，我們先考慮將最後一個拿下來（沒辦法，奇葩的規則，只有拿下前n-1個且保留第n個，才能拿下第n+1個）。那麼先將前n-1個取下來，需要f(n-1)步；再將第n+1個取下來，需要1步；再把前n-1個放回去，需要f(n-1)步（這裡需要解釋一下，前n-1個必須放回去，不然第n個無法取下，而且取下需要f(n-1)步，那麼放回去也需要f(n-1)，直接相反的步驟再來一遍，不明白的多讀幾遍題目）；最後，將前n個都取下，需要f(n)步。所以得到遞推式 f(n+1)=f(n)+2f(n-1)+1。對於遞推式中有常數項，只需要給出一行或列留給常數項。
構造矩陣
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 200907;//1e9+7;
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix& rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
int n;
matrix ant;
ant.init(3);
ant.maze[0][0]=ant.maze[0][1]=ant.maze[2][0]=ant.maze[2][2]=1;
ant.maze[1][0]=2;
while(scanf("%d",&n),n)
{
if(n<=2)
{
if(n==1)printf("%d\n",1);
else printf("%d\n",2);
continue;
}
matrix ans;
ans.init(3);
ans.maze[0][0]=2;
ans.maze[0][1]=ans.maze[0][2]=1;
ans=ans*qlow(ant,n-2);
printf("%I64d\n",ans.maze[0][0]);
}
return 0;
}第七題 hdu-2276
分析：對於01串，如果前一位的數位為1，那麼這一位就翻轉（0變1，1變0）；前一位為0，那麼這一位就不變。這就像前一位與這一位相加，並且mod2，如果前一位是0，那麼這一位加上前一位再mod2數值不會改變；如果前一位是1，那麼這一位就相當於+1mod2，正好和翻轉的定義相同。這樣我們就可以構造出矩陣
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 100+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;
char str[maxn];
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix & rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%2;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
int m;
while(~scanf("%d",&m))
{
scanf("%s",str);
int n=strlen(str);
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
if(str[i]=='1')ans.maze[0][i]=1;
matrix ant;
ant.init(n);
for(int i=0;i<n;i++)
{
ant.maze[i][i]=1;
if(i==0)ant.maze[n-1][0]=1;
else ant.maze[i-1][i]=1;
}
ant=qlow(ant,m);
ans=ans*ant;
for(int i=0;i<n;i++)
printf("%d",ans.maze[0][i]);
puts("");
}
return 0;
}第八題 hdu-2254
分析：學過離散數學的都知道，如果用鄰接矩陣來存圖的話，矩陣A的n次方里面的（A^n）[i][j]表示從i到j且路徑長度為n的不同路徑條數。這一定理正好適用這一題目。我用鄰接矩陣來建圖（單向邊，一條邊對應一個++操作）得到矩陣A，題目要求的就是｛A^(t1)｝[v1][v2]+｛A^(t1+1)｝[v1][v2]+｛A^(t1+2)｝[v1][v2]+......+｛A^(t2)｝[v1][v2]的和取模。但是如果一個一個算的話特別麻煩，還有超時的可能，所以我們進一步改進思路。
對於，我們不妨設Z(n)=A^1+A^2+...+A^n
那麼可以構造遞推關係 A^(n+1)=A^n * A Z(n+1)=Z(n)+ A^(n+1) =Z(n)+A^n * A
構造矩陣
但是這道題目略坑，有幾個注意點：
1、城市的編號不是從0到n-1，而是隨便的一個數字，需要離散化否則不能存相關資訊
2、城市數不超過30，也就是說我的方法開矩陣不超過60，但是我殘念的一開始以為最多可能有20000個不同城市 血崩！
3、圖中可能有重邊，所以別用=1，要用++操作
4、詢問中v1，v2可能在前面的城市編號集中沒有出現，那麼此時答案為0
5、t1可能比t2大，這種情況你就交換下t1，t2
這道題要做對也是真心不容易！
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 60+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 2008;//1e9+7;
ll gra[10000+5][2];
ll num[20000+5];
int sum;
int bs(ll k)
{
int m,l=0,r=sum-1;
while(l<=r)
{
m=(l+r)>>1;
if(k==num[m])return m;
else if(k>num[m])l=m+1;
else r=m-1;
}
return -1;
}
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix & rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
int n;
while(~scanf("%d",&n))
{
sum=n*2;
for(int i=0;i<n;i++)
{
scanf("%lld %lld",&gra[i][0],&gra[i][1]);
num[i<<1]=gra[i][0];
num[i<<1|1]=gra[i][1];
}
sort(num,num+sum);
sum=unique(num,num+sum)-num;
matrix ans;
ans.init(sum<<1);
matrix ant;
ant.init(sum<<1);
int a,b;
for(int i=0;i<n;i++)
{
a=bs(gra[i][0]);
b=bs(gra[i][1]);
ans.maze[a][b]++;
ans.maze[a][b+sum]++;
}
for(int i=0;i<sum;i++)
ans.maze[i+sum][i+sum]=1;
int q;
scanf("%d",&q);
while(q--)
{
int t1,t2;
ll v1,v2;
scanf("%lld %lld %d %d",&v1,&v2,&t1,&t2);
v1=bs(v1);
v2=bs(v2);
if(v1==-1||v2==-1)
{
puts("0");
continue;
}
if(t1>t2)swap(t1,t2);
t1--;
int d,d1,d2;
if(t1<=0)d1=0;
else
{
ant=qlow(ans,t1);
d1=ant.maze[v1][v2+sum];
}
if(t2<=0)d2=0;
else
{
ant=qlow(ans,t2);
d2=ant.maze[v1][v2+sum];
}
d=d2-d1;
while(d<0)d+=mod;
printf("%d\n",d);
}
}
return 0;
}第九題 hdu-3117
分析：其實這道題不應該放在矩陣快速冪專題的，因為這道題關於矩陣的部分很簡單-。-，難的是題目的數學部分-。-
求解斐波那契數列的值，對於數值位不超過8位的直接輸出數值；但是對於大於8位的，輸出值的前四位和後四位（格式自己看-。-）。後四位當然非常簡單了，f(n)=f(n-1)+f(n-2) ,mod 10000,我就不再贅述了
但是對於前四位還是挺麻煩的，這裡要用到對數的優良性質。logA(b^c)=c*logA(c) logA(b*c)=logA(b)+logA(c)
對於一個大數（不如就選10234432）,怎麼知道它的前四位呢
log10(10234432)=log10(1.0234432*10^7)=log10(1.0234432)+7
7是log10(10234432)的整數部分，log10(1.0234432)是log10(10234432)的小數部分
可以明顯看出pow(10,log10(1.0234432))=1.0234432 從而可以簡單地得到了前幾個位數
有人會問為什麼要這麼麻煩，不如以字串的形式讀入然後輸出前四位不就好了嗎 ？
對！沒錯，對於確定的大數，確實是可以的，但是對b^c來說卻不能使用，但是取對數的方法卻可以
對於這道題 我們考慮斐波那契數列的通項（有人問，比賽的時候你也知道通項？我的回答是可以，給我類似斐波那契的遞推式我就可以求出其通項，這個在數學競賽中有學到，請讀者自行百度）
F(n)=(1/sqrt(5))*{[(1+sqrt(5))/2]^n-[(1-sqrt(5))/2]^n}=(1/sqrt(5))*{[(1+sqrt(5))/2]^n}*{1-[(1-sqrt(5))/(1+sqrt(5))]^n}
設p=(1-sqrt(5))/2
s=log10(F(n))=-1/2*log10(5)+n*log10(p)+Δ（Δ非常小，可以忽略不計，Δ小於log10(1)）
s去掉它的整數部分，然後(int) [ pow（10，s-（int）s）*1000 ] 就得到了結果
我只是簡單搬運而且=、=
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 60+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 10000;//1e9+7;
int fa[60];
struct matrix
{
int n;
ll maze[5][5];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix & rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
return ans;
}
};
matrix qlow(matrix a,int n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
fa[0]=0;
fa[1]=1;
for(int i=2;i<40;i++)
fa[i]=fa[i-1]+fa[i-2];
int n;
while(~scanf("%d",&n))
{
if(n<40)
{
printf("%d\n",fa[n]);
continue;
}
double p=(sqrt(5.0)+1.0)/2.0;
double s=1.0*n*log10(p)-log10(5.0)/2;
s=s-(int)s;
int anw=(int)(pow(10.0,s)*1000);
printf("%d...",anw);
matrix ant;
ant.init(2);
ant.maze[0][0]=ant.maze[0][1]=ant.maze[1][0]=1;
ant=qlow(ant,n-1);
printf("%04d\n",ant.maze[0][0]);
}
return 0;
}第十題 hdu-4686
分析：這是一道顯然的矩陣題（刷完前面那麼多道，如果還沒有一眼看出來的話，建議停下來重新理解演算法-。-）
感覺難度不是很大，但是這道題有多個遞推關係，而且還涉及到常數項和兩個不相關的遞推關係乘積，稍微有點麻煩。但是冷靜下來，對每一個遞推關係進行分析和對你所要求的項進行分解。
由AoD(n)=a(0)b(0)+a(1)b(1)+...+a(n-1)b(n-1) 可得 AoD(n+1)=AoD(n)+a(n)b(n)
那麼a(n)b(n)怎麼求呢 現在只知道 a(n)=a(n-1)*ax+ay b(n)=b(n-1)*bx+by
那麼可不可以利用已經得到的a(n-1)和b(n-1)和a(n-1)b(n-1)
可以得到 a(n)b(n)=ax*bx*a(n-1)*b(n-1)+ax*by*a(n-1)+ay*bx*b(n-1)+ay*by
那麼問題又來了 a(n)和 b(n)又怎麼單獨求呢 這個題目已經給出 a(n)=a(n-1)*ax+ay b(n)=b(n-1)*bx+by
分析清楚了所有的遞推關係後，可以來構造矩陣了
注意點：a0，ax，ay，b0，bx，by都可能大於1e7+9，所以需要先取一次模
#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 1+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;
ll mul(ll a,ll b)
{
ll ans=0;
while(b)
{
if(b&1)ans=(ans+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return ans;
}
struct matrix
{
int n;
ll maze[maxn][maxn];
void init(int n)
{
this->n=n;
clr(maze,0);
}
matrix operator * (const matrix & rhs)
{
matrix ans;
ans.init(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
return ans;
}
};
matrix qlow(matrix a,ll n)
{
matrix ans;
ans.init(a.n);
for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
while(n)
{
if(n&1)ans=ans*a;
a=a*a;
n>>=1;
}
return ans;
}
int main()
{
//freopen("d:\\acm\\in.in","r",stdin);
ll n;
while(~scanf("%lld",&n))
{
ll a0,b0,ax,bx,ay,by;
scanf("%lld %lld %lld %lld %lld %lld",&a0,&ax,&ay,&b0,&bx,&by);
a0%=mod;
b0%=mod;
ax%=mod;
bx%=mod;
ay%=mod;
by%=mod;
if(n==0)
{
puts("0");
continue;
}
matrix ans;
ans.init(5);
ans.maze[0][0]=1;
ans.maze[0][1]=a0;
ans.maze[0][2]=b0;
ans.maze[0][3]=a0*b0%mod;
matrix ant;
ant.init(5);
ant.maze[0][0]=1;
ant.maze[0][1]=ay;
ant.maze[1][1]=ax;
ant.maze[0][2]=by;
ant.maze[2][2]=bx;
ant.maze[0][3]=ay*by%mod;
ant.maze[1][3]=ax*by%mod;
ant.maze[2][3]=ay*bx%mod;
ant.maze[3][3]=ax*bx%mod;
ant.maze[3][4]=ant.maze[4][4]=1;
ant=qlow(ant,n);
ans=ans*ant;
printf("%lld\n",ans.maze[0][4]);
}
return 0;
}這裡先小結一下，矩陣快速冪型別的題目非常像科學歸納法（第一科學歸納法、第二科學歸納法還有其他型別的），主要是求得當前項與前幾項或者常數項的遞推關係，從而構造矩陣求解。