矩陣快速冪專題(一)
最近閒來無事,準備集中精力刷一波數論與圖論。矩陣快速冪是數論裡面的重要組成部分,值得我好好學習一下。因為題目比較多,分析也比較多,所以將此專題分成幾個部分。做完這一專題,可能會暫時轉向圖論部分,然後等我組合數學學得差不多了,再回過頭來繼續做數論題。
矩陣快速冪演算法的核心思想是將問題建模轉化為數學模型(有一些簡單題目是裸的矩陣模型,但是大部分難題就是難在要構造矩陣,用矩陣方法解決問題),推倒遞推式,構造計算矩陣,用快速冪的思想求解矩陣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;
}
第三題 hdu-2604
分析:我這裡只提出我的方法,其他網上大神的方法請到別的部落格上查詢。
對於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;
}
第六題 hdu-2842分析:題目講的好像是九連環(好像又不是,英語不好),-。- 可憐我並沒有玩過
這道題非常像漢諾塔,但是題目只是求最小步數,而不需要求出具體的步驟。我假設對於長度為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;
}
這裡先小結一下,矩陣快速冪型別的題目非常像科學歸納法(第一科學歸納法、第二科學歸納法還有其他型別的),主要是求得當前項與前幾項或者常數項的遞推關係,從而構造矩陣求解。