写的不错的博客：http://www.cnblogs.com/yan-boy/archive/2012/11/29/2795294.html
优点：根据数列递推式快速计算数列an的值（当n很大时），an可以是数，也可以是矩阵
步骤：由数列递推式构造矩阵，然后用矩阵快速幂计算矩阵的幂。
构造矩阵：
对于an =x*an-1 +y*an-2 ,可以构造矩阵为：
[an ,an-1]=[an-1 ,an-2]*A
A=[x 1
y 0];
矩阵快速幂模板：
（2018.8.14修改）
struct Matrix { int a[3][3]; void init() { for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) a[i][j] = 0; } } void _init() { init(); for (int i = 0; i < 3; i++) a[i][i] = 1; } }A, B; Matrix mul(Matrix a, Matrix b) { Matrix ans; ans.init(); for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { if(a.a[i][j]) { for (int k = 0; k < 3; k++) ans.a[i][k] = (ans.a[i][k] + 1LL * a.a[i][j] * b.a[j][k]) % MOD; } } } return ans; } Matrix q_pow(Matrix a, int k) { Matrix ans; ans._init(); if(k <= 0) return ans; while(k) { if(k&1) ans = mul(ans, a); a = mul(a, a); k >>= 1; } return ans; }
代码：
#include<iostream> #include<cstring> #include<cstdio> #include<algorithm> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) const int MOD=10000; struct Matrix { ll a[2][2]; }A,res; Matrix multiply(Matrix x,Matrix y) { Matrix res; mem(res.a,0); 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 qpow(int n) { while(n) { if(n&1)res=multiply(res,A); A=multiply(A,A); n=n>>1; } } int main() { ios::sync_with_stdio(false); cin.tie(0); int n; while(cin>>n) { if(n==-1)break; mem(res.a,0); for(int i=0;i<2;i++)res.a[i][i]=1; A.a[0][0]=1; A.a[0][1]=1; A.a[1][0]=1; A.a[1][1]=0; qpow(n); cout<<res.a[0][1]<<endl; } return 0; }
代码：
#include<bits/stdc++.h> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) const int MOD=998244353; struct Matrix { ll a[2][2]; }res,A; Matrix mul(Matrix x,Matrix y) { Matrix ans; mem(ans.a,0); for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++) ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%MOD; return ans; } void qpow(int n) { mem(res.a,0); for(int i=0;i<2;i++) res.a[i][i]=1; A.a[0][0]=A.a[0][1]=A.a[1][0]=1; A.a[1][1]=0; while(n) { if(n&1)res=mul(res,A); A=mul(A,A); n>>=1; } } int main() { ios::sync_with_stdio(false); cin.tie(0); int k; while(cin>>k) { int n=k*2+3; qpow(n); cout<<res.a[0][1]-1<<endl; } return 0; }
方法1：
构造矩阵：
B=(A E
O E)
Bk+1=(Ak+1 E+A+A2 +A3 +...+Ak
O E )
代码1：
#include<iostream> #include<cstring> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) int m,n; struct Matrix { ll a[65][65]; }res,A; Matrix mul(Matrix x,Matrix y) { Matrix ans; mem(ans.a,0); for(int i=0;i<2*n;i++) for(int j=0;j<2*n;j++) for(int k=0;k<2*n;k++) ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%m; return ans; } void qpow(int k) { mem(res.a,0); for(int i=0;i<2*n;i++)res.a[i][i]=1; for(int i=0;i<n;i++)A.a[i][i+n]=1; for(int i=n;i<2*n;i++)A.a[i][i]=1; while(k) { if(k&1)res=mul(res,A); A=mul(A,A); k>>=1; } } int main() { ios::sync_with_stdio(false); cin.tie(0); int k; cin>>n>>k>>m; for(int i=0;i<n;i++) for(int j=0;j<n;j++) cin>>A.a[i][j]; qpow(k+1); for(int i=0;i<n;i++) { for(int j=n;j<2*n;j++) { if(j!=n+i)cout<<res.a[i][j]; else cout<<(res.a[i][j]-1+m)%m; if(j!=2*n-1)cout<<' '; } cout<<endl; } return 0; }
方法2：
由Fn =A*Fn-1 +A构造矩阵：
(Fn ,E)=(Fn-1,E)*B
B=(A O
A E)
(Fn ,E)=(F1,E)*Bn-1
代码2：
#include<iostream> #include<cstring> #include<cstdio> #include<algorithm> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) int n,m,k; struct Matrix { ll a[65][65]; }A,res,t; Matrix mul(Matrix x,Matrix y) { Matrix ans; mem(ans.a,0); for(int i=0;i<2*n;i++) for(int j=0;j<2*n;j++) for(int k=0;k<2*n;k++) ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%m; return ans; } void qpow(int k) { mem(res.a,0); mem(A.a,0); for(int i=0;i<2*n;i++)res.a[i][i]=1; for(int i=0;i<n;i++) for(int j=0;j<n;j++) A.a[i][j]=A.a[i+n][j]=t.a[i][j]; for(int j=n;j<2*n;j++)A.a[j][j]=1; while(k) { if(k&1)res=mul(res,A); A=mul(A,A); k=k>>1; } } int main() { ios::sync_with_stdio(false); cin.tie(0); cin>>n>>k>>m; for(int i=0;i<n;i++) for(int j=0;j<n;j++) cin>>t.a[i][j]; qpow(k-1); for(int i=0;i<n;i++) t.a[i][i+n]=1; /*for(int i=0;i<2*n;i++) for(int j=0;j<2*n;j++) { cout<<res.a[i][j]<<' '; if(j==2*n-1)cout<<endl; }*/ for(int i=0;i<n;i++) { for(int j=0;j<n;j++) { int ans=0; for(int k=0;k<2*n;k++) { ans=(ans+t.a[i][k]*res.a[k][j])%m; } cout<<ans; if(j!=n-1)cout<<' '; } cout<<endl; } return 0; }
写的不错的博客：http://blog.csdn.net/magicnumber/article/details/6217602
要用到一个对于稀疏矩阵乘法的加速技巧
代码:
#include<iostream> #include<cstring> #include<cstdio> #include<algorithm> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) int n,m,k; struct Matrix { ll a[105][105]; }A,res,unit; Matrix mul(Matrix x,Matrix y) { Matrix ans; mem(ans.a,0); for(int i=0;i<n+1;i++) for(int j=0;j<n+1;j++) { if(x.a[i][j]) for(int k=0;k<n+1;k++) ans.a[i][k]+=x.a[i][j]*y.a[j][k]; /*for(int k=0;k<n+1;k++) ans.a[i][j]+=x.a[i][k]*y.a[k][j];*/ } return ans; } void init() { mem(unit.a,0); for(int i=0;i<n+1;i++) { unit.a[i][i]=1; } } void qpow(int k) { init(); res=unit; while(k) { if(k&1)res=mul(res,A); A=mul(A,A); k>>=1; } } int main() { ios::sync_with_stdio(false); cin.tie(0); while(cin>>n>>m>>k) { if(n==0&&m==0&&k==0)break; int i,j; mem(A.a,0); char s; init(); res=unit; while(k--) { cin>>s; if(s=='g') { cin>>i; init(); unit.a[i-1][n]=1; res=mul(unit,res); } else if(s=='e') { cin>>i; init(); unit.a[i-1][i-1]=0; res=mul(unit,res); } else if(s=='s') { cin>>i>>j; init(); unit.a[i-1][i-1]=0; unit.a[i-1][j-1]=1; unit.a[j-1][j-1]=0; unit.a[j-1][i-1]=1; res=mul(unit,res); } } A=res; qpow(m); for(int i=0;i<n;i++) { cout<<res.a[i][n]; if(i!=n-1)cout<<' '; else cout<<endl; } } return 0; }
例题5：HDU - 5015
首先，如果没有最上面一行，每新增加一列都是求一遍前缀和，那么矩阵就是一个下三角全为1的矩阵
但是，这道题最上面有数，那我们矩阵稍微改一改
10 0 0 0 ... 1 23
10 1 0 0 ... 1 a1
10 1 1 0 ... 1 a2
10 1 1 1 ... 1 * a3
...... .
10 1 1 1 ... 1 an
0 0 0 0 ... 1 3
这个相当于算第一列，算m列就是乘m次左边的那个矩阵
#pragma GCC optimize(2) #pragma GCC optimize(3) #pragma GCC optimize(4) #include<bits/stdc++.h> using namespace std; #define fi first #define se second #define pi acos(-1.0) #define LL long long //#define mp make_pair #define pb push_back #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r #define ULL unsigned LL #define pll pair<LL, LL> #define pii pair<int, int> #define piii pair<pii, int> #define mem(a, b) memset(a, b, sizeof(a)) #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); #define fopen freopen("in.txt", "r", stdin);freopen("out.txt", "w", stout); //head const int MOD = 1e7 + 7; int a[100], n; struct Matrix { int a[15][15]; void init() { for (int i = 0; i <= n+1; i++) { for (int j = 0; j <= n+1; j++) { a[i][j] = 0; } } } void _init() { init(); for (int i = 0; i <= n+1; i++) a[i][i] = 1; } }T; Matrix mul(Matrix A, Matrix B) { Matrix ans; ans.init(); for (int i = 0; i <= n+1; i++) { for (int j = 0; j <= n+1; j++) { if(A.a[i][j]) for (int k = 0; k <= n+1; k++) ans.a[i][k] = (ans.a[i][k] + 1LL * A.a[i][j] * B.a[j][k]) % MOD; } } return ans; } Matrix q_pow(Matrix A, int k) { Matrix ans; ans._init(); while(k) { if(k&1) ans = mul(ans, A); A = mul(A, A); k >>= 1; } return ans; } void init() { T.a[0][0] = 10; for (int i = 1; i <= n; i++) T.a[i][0] = 10; for (int i = 1; i <= n; i++) { for (int j = 1; j <= i; j++) T.a[i][j] = 1; } for (int i = 0; i <= n+1; i++) T.a[i][n+1] = 1; } int main() { int m; while(~ scanf("%d %d", &n, &m)) { for (int i = 1; i <= n; i++) scanf("%d", &a[i]); T.init(); init(); T = q_pow(T, m); int ans = T.a[n][0] * 23LL % MOD; for (int i = 1; i <= n; i++) ans = (ans + 1LL*T.a[n][i]*a[i]) % MOD; ans = (ans + 1LL*T.a[n][n+1]*3LL) % MOD; printf("%d ", ans); } return 0; }