矩阵乘法是一种高效的算法可以把一些一维递推优化到log（ n )，还可以求路径方案等，所以更是是一种应用性极强的算法。矩阵，是线性代数中的基本概念之一。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。由于它把许多数据紧凑的集中到了一起，所以有时候可以简便地表示一些复杂的模型。矩阵乘法看起来很奇怪，但实际上非常有用，应用也十分广泛。
基本定义
它是这样定义的，只有当矩阵A的列数与矩阵B的行数相等时A×B才有意义。一个m×n的矩阵a（m,n)左乘一个n×p的矩阵b（n,p)，会得到一个m×p的矩阵c（m,p)，满足
矩阵乘法满足结合率，但不满足交换率
一般的矩乘要结合快速幂才有效果``
http://poj.org/problem?id=3070
#include<cstdio> #include<iostream> using namespace std; const int MOD = 10000; struct matrix{ int m[2][2]; }ans, base; matrix mul(matrix a, matrix b) { matrix tmp; for(int i=0; i<2; ++i) { for(int j=0; j<2; ++j) { tmp.m[i][j] = 0; for(int k=0; k<2; ++k) tmp.m[i][j] = (tmp.m[i][j]+a.m[i][k]*b.m[k][j])%MOD; } } return tmp; } int fast_mod(int n) { base.m[0][0] = base.m[0][1] = base.m[1][0] = 1; base.m[1][1] = 0; ans.m[0][0] = ans.m[1][1] = 1; ans.m[0][1] = ans.m[1][0] = 0; while(n) { if(n&1) { ans = mul(ans, base); } base = mul(base, base); n >>= 1; } return ans.m[0][1]; } int main() { int n; while(scanf("%d", &n)&&n!=-1) { printf("%d ", fast_mod(n)); } return 0; }
稍加优化：
#include<cstdio> #include<cstring> #include<iostream> using namespace std; const int MOD = 10000; struct matrix{ int m[2][2]; }ans, base; matrix mul(matrix a, matrix b) { matrix tmp; memset(tmp.m, 0, sizeof(tmp.m)); for(int i=0; i<2; i++) for(int k=0; k<2; k++) if(a.m[i][k]) for(int j=0; j<2; j++) tmp.m[i][j] = (tmp.m[i][j]+a.m[i][k]*b.m[k][j])%MOD; return tmp; } int fast_mod(int n) { base.m[0][0] = base.m[0][1] = base.m[1][0] = 1; base.m[1][1] = 0; ans.m[0][0] = ans.m[1][1] = 1; ans.m[0][1] = ans.m[1][0] = 0; while(n) { if(n&1) { ans = mul(ans, base); } base = mul(base, base); n >>= 1; } return ans.m[0][1]; } int main() { int n; while(scanf("%d", &n)&&n!=-1) { printf("%d ", fast_mod(n)); } return 0; }
《2》http://www.lightoj.com/volume_showproblem.php?problem=1096
自行构造矩阵吧！
#include<iostream> #include<cstdio> #include<cstring> using namespace std; int k,M=10007; struct mat{ int ma[4][4]; mat(){memset(ma, 0, sizeof(ma));} mat(int a){ memset(ma, 0, sizeof(ma)); for(int i=0; i<4; i++) ma[i][i] = 1; } void operator *= (const mat &x) { mat tmp; for(k=0; k<4; k++) for(int i=0; i<4; i++) if(ma[i][k]) for(int j=0; j<4; j++) tmp.ma[i][j] = (tmp.ma[i][j]+ma[i][k]*x.ma[k][j])%M; for(int i=0; i<4;i++) for(int j=0; j<4; j++) ma[i][j] = tmp.ma[i][j]; } }; mat pow_mod(mat a, int b) { mat tmp = mat(1); while(b>0) { if(b&1) tmp*=a; b>>=1; a*=a; } return tmp; } int main() { int a, b, c, n, T; int ans = 0; scanf("%d", &T); for(int kk=1; kk<=T; kk++) { scanf("%d%d%d%d", &n, &a, &b, &c); mat A; A.ma[0][0] = a; A.ma[0][2] = b; A.ma[0][3] = 1; A.ma[1][0] = 1; A.ma[2][1] = 1; A.ma[3][3] = 1; if(n>3) { A = pow_mod(A, n-3); ans = (A.ma[0][0]*c+A.ma[0][3]*c)%M; } else if(n==3) ans = c%M; else ans = 0; printf("Case %d: %d ", kk, ans); } return 0; }
《3》来个较难的矩阵快速幂。
POJ3233
http://poj.org/problem?id=3233
题目大意：给定矩阵A，求A + A^2 + A^3 + ... + A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据mod m。k<=10^9。
这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：
A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)
应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A + A^2 + A^3，即可得到原问题的答案。
k=7时，有：A+A^2+A^3+A^4+A^5+A^6+A^7 =A+A^2+A^3+A^3*(A+A^2+A^3+A^3*A^1)
进行了两次二分思想， 绝对是很巧妙地， 代码实现的也很精湛。你们不用怀疑这是草滩小恪在自卖自夸， 因为草滩小恪只想到了两次二分的思路， 并不能用代码实现。 但是草滩小恪在百小度的帮忙下看到的一个巨巨的博客http://www.cnblogs.com/DreamUp/archive/2010/07/27/1786225.html。
草滩小恪“偷”来的代码：
#include<stdio.h> #include<string.h> int n,kk,m; int ans[31][31],a[31][31],ori[31][31]; int tmp[31][31],tp[31][31]; void calSqu()//矩阵 A = A*A。 { int i,j,k; memset(tmp,0,sizeof(tmp)); for(i=0;i<n;i++) for(k=0;k<n;k++) if(a[i][k]) for(j=0; j<n; j++) tmp[i][j]+=(a[i][k]*a[k][j])%m; for(i=0;i<n;i++) for(j=0;j<n;j++) a[i][j]=tmp[i][j]%m; } void cal(int t) { int i,j,k; if(t==1) return; if(t%2==0) { cal(t/2); memset(tmp,0,sizeof(tmp)); for(i=0;i<n;i++) for(k=0;k<n;k++) if(a[i][k]) for(j=0;j<n;j++) tmp[i][j]+=(a[i][k]*ans[k][j])%m; for(i=0;i<n;i++) for(j=0;j<n;j++) ans[i][j]=(ans[i][j]+tmp[i][j])%m; calSqu(); } else{ cal(t/2); memset(tp,0,sizeof(tp)); for(i=0;i<n;i++) for(j=0;j<n;j++) tp[i][j]=ans[i][j]; for(i=0;i<n;i++) for(j=0;j<n;j++) { for(k=0;k<n;k++) tp[i][j]+=(a[i][k]*ori[k][j])%m; tp[i][j]%=m; } memset(tmp,0,sizeof(tmp)); for(i=0;i<n;i++) for(j=0;j<n;j++) { for(k=0;k<n;k++) tmp[i][j]+=(a[i][k]*tp[k][j])%m; tmp[i][j]%=m; } for(i=0;i<n;i++) for(j=0;j<n;j++) ans[i][j]=(ans[i][j]+tmp[i][j])%m; calSqu(); memset(tmp,0,sizeof(tmp)); for(i=0;i<n;i++) for(j=0;j<n;j++) { for(k=0;k<n;k++) tmp[i][j]+=(a[i][k]*ori[k][j])%m; tmp[i][j]%=m; } for(i=0;i<n;i++) for(j=0;j<n;j++) a[i][j]=tmp[i][j]; } } int main() { int i,j; scanf("%d%d%d",&n,&kk,&m); for(i=0;i<n;i++) for(j=0;j<n;j++) { scanf("%d",&a[i][j]); ans[i][j]=ori[i][j]=a[i][j]; } cal(kk); for(i=0;i<n;i++) { for(j=0;j<n-1;j++) printf("%d ",ans[i][j]%m); printf("%d ",ans[i][j]%m); } return 0; }
不过上面的那个代码有些难懂， 不是太清晰。
思路：
(A + A2 + … + Ak/2) + Ak/2 *(A + A2 + … + Ak/2)
S(k) =
(A + A2 + … + Ak/2) + Ak/2 *(A + A2 + … + Ak/2) + Ak
即：
S(k/2) + Ak/2* S(k/2) n为偶数
S(k) =
S(k/2) + Ak/2* S(k/2) + Ak n为奇数
来个更优的解法。
构造一个分块矩阵（矩阵的元素也为矩阵）
设Sk = I + A + ……+ Ak-1
则有Sk+1 = Sk + Ak
=>
|Sk+1| | 1 1 | |Sk|
|Ak+1| = | 0 A | |Ak|
这种方法巧妙地应用的矩阵的乘法， 有矩阵乘法的结合律， 就把问题转化成了 --------> 矩阵快速幂。代码实现中有细节问题， 请细加揣摩。（不妨在练习本上模拟一下）。
#include <cstdio> #include <cstdlib> #include <sstream> #include <iostream> #include <cmath> #include <cstring> #include <algorithm> #include <string> #include <utility> #include <vector> #include <queue> #include <map> #include <set> using namespace std; typedef long long ll; typedef pair<int,int> PII; int n,k,M; struct mat{ int ma[62][62]; mat(){memset(ma,0,sizeof(ma));}//¹¹ÔìÁã¾ØÕó mat(int a){//¹¹Ô쵥λ¾ØÕó memset(ma,0,sizeof(ma)); for(int i=0; i<n; i++) ma[i][i] = 1; } void operator *= (const mat &x) {//¾ØÕó³Ë·¨ mat tmp; for(k=0; k<n; k++) for(int i=0; i<n; i++) if(ma[i][k]){ for(int j=0; j<n; j++) { tmp.ma[i][j] += ma[i][k]*x.ma[k][j]; tmp.ma[i][j] %= M; } } for(int i=0; i<n; i++) for(int j=0; j<n; j++) ma[i][j] = tmp.ma[i][j]; } }; mat pow_mod(mat a,int b){//¾ØÕó¿ìËÙÃÝ mat tmp = mat(1); while(b>0){ if(b&1)tmp*=a; b >>= 1; a*=a; } return tmp; } int main(){ while(~scanf("%d%d%d",&n,&k,&M)){ mat A; for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { scanf("%d",&A.ma[i+n][j+n]); } A.ma[i][i] = A.ma[i][i+n] = 1; } n <<= 1; mat ans = pow_mod(A,k+1); n >>= 1; for(int i=0; i<n; i++) { for(int j=0; j<n; j++) { if(i==j)ans.ma[i][j+n] = (ans.ma[i][j+n]-1+M)%M; printf("%d%c",ans.ma[i][j+n],(j==n-1)?' ':' '); } } } return 0; }
《4》http://www.lightoj.com/volume_showproblem.php?problem=1244
这道题的递推式的发现不是那么明显哦（嘿嘿！）。
设f(n)为2*n铺满的方案数。 考虑最后一块砖的放法，有四种情况：
其中两种可以由f(n-2)和f(n-1)得来，另外两种则不能，所以需要加多两种状态。
设u(n)为放了n列但右上角空着的方案数；
v(n)为放了n列但右下角空着的方案数。
可以得到f(n) = f(n-1) + f(n-2) + u(n-1) + v(n-1)
对于u(n)，同样考虑最后一块的放法：
则u(n) = v(n-1) + f(n-2)
同理有v(n) = u(n-1) + f(n-2)
最后构造矩阵，矩阵快速幂。
只要推出递推式， 矩阵就很好构造啦！。 不得不说草滩小恪这次还真是机智！，----难得机智一回。
#include<iostream> #include<cstdio> #include<cstring> using namespace std; const int M = 10007; struct mat{ int ma[4][4]; mat(){memset(ma, 0, sizeof(ma));} mat(int a){ memset(ma, 0, sizeof(ma)); for(int i=0; i<4; i++) ma[i][i] = 1; } void operator *= (const mat &x) { mat tmp; for(int k=0; k<4; k++) for(int i=0; i<4; i++) if(ma[i][k]) for(int j=0; j<4; j++) tmp.ma[i][j] = (tmp.ma[i][j]+ma[i][k]*x.ma[k][j])%M; for(int i=0; i<4;i++) for(int j=0; j<4; j++) ma[i][j] = tmp.ma[i][j]; } }; mat pow_mod(mat a, int b) { mat tmp = mat(1); while(b>0) { if(b&1) tmp*=a; b>>=1; a*=a; } return tmp; } int main() { int n, T; int ans = 0; scanf("%d", &T); for(int kk=1; kk<=T; kk++) { scanf("%d", &n); mat A; A.ma[0][0] = 1; A.ma[0][1] = 1; A.ma[0][2] = 1; A.ma[0][3] = 1; A.ma[1][0] = 1; A.ma[2][1] = 1; A.ma[2][3] = 1; A.ma[3][1] = 1; A.ma[3][2] = 1; if(n>2) { A = pow_mod(A, n-2); ans=(2*A.ma[0][0]+A.ma[0][1]+A.ma[0][2]+A.ma[0][3])%M; } else ans = n; printf("Case %d: %d ", kk, ans); } return 0; }
总结： 以上问题其实都是----------------->>>矩阵解常系数线性递推方程
递推方程Fn = aFn-1 + bFn-2 + cFn-3 …。。。。。。。。。。。。 可以写成两个向量相乘
(Fn-1 , Fn-2 , Fn-3 , … )*(a , b , c , … )然而十分巧合的事情来啦（哈哈！）。细心观察你就会惊奇的发现这个式子和跟矩阵乘法的C[i,j] = ∑A[i,k]*B[k,j]很是神似。
于是乎可以对向量T=(Fn-1 , Fn-2 , Fn-3 , … )构造一个变换矩阵A ，使得T’= T*A
然后进行迭代， 然后你就会又惊奇的发现， A*A*A，，，，，重复乘了好多的A哦！ 。 毫无疑问， 这时你想到了矩阵快速幂， 于是你就把这个问题解决啦。 不用谢小恪， 小恪是雷锋！
继续练兵：《5》http://acm.hdu.edu.cn/showproblem.php?pid=5171
看到题分析了一下， 直接暴力了一下， 果断的超时啦！ 暴力代码如下：
#include<iostream> using namespace std; const int maxn = 100000+5; const int MOD = 10000007; int a[maxn]; int main() { int max1, max2, flag, wap; int n, k; int ans=0; while(scanf("%d%d", &n, &k)==2) { max1=0; max2 = 0; flag=0; ans = 0; for(int i=0; i<n; i++) { scanf("%d", &a[i]); if(a[i]>max2) { max2 = a[i]; flag=i; } ans=(ans+a[i])%MOD; } for(int i=0; i<n; i++) if(a[i]>max1&&i!=flag) { max1 = a[i]; } for(int i=0; i<k; i++) { ans=(ans+max1+max2)%MOD; wap = (max1+max2)%MOD; max1 = max2; max2 = wap; } printf("%d ", ans); } return 0; }
然后构造了一个矩阵快速幂算法， 结果还是超时啦！！！（哭！）
#include<iostream> #include<cstring> using namespace std; const int maxn = 100000+5; const int MOD = 10000007; int a[maxn]; struct mat{ int ma[3][3]; mat(){memset(ma,0,sizeof(ma));} mat(int a){ memset(ma,0,sizeof(ma)); for(int i=0; i<3; i++) ma[i][i] = 1; } void operator *= (const mat &x) {//¾ØÕó³Ë·¨ mat tmp; for(int k=0; k<3; k++) for(int i=0; i<3; i++) if(ma[i][k]){ for(int j=0; j<3; j++) { tmp.ma[i][j] += ma[i][k]*x.ma[k][j]; tmp.ma[i][j] %= MOD; } } for(int i=0; i<3; i++) for(int j=0; j<3; j++) ma[i][j] = tmp.ma[i][j]; } }; mat pow_mod(mat a,int b){//¾ØÕó¿ìËÙÃÝ mat tmp = mat(1); while(b>0){ if(b&1)tmp*=a; b >>= 1; a*=a; } return tmp; } int main() { long long max1, max2, flag, wap; long long n, k; long long ans=0; while(scanf("%d%d", &n, &k)==2) { max1=0; max2 = 0; flag=0; ans = 0; for(int i=0; i<n; i++) { scanf("%d", &a[i]); if(a[i]>max2) { max2 = a[i]; flag=i; } ans=(ans+a[i])%MOD; } for(int i=0; i<n; i++) if(a[i]>max1&&i!=flag) { max1 = a[i]; } mat A; A.ma[0][0] = 1; A.ma[0][1] = 1; A.ma[0][2] = 1; A.ma[1][1] = 1; A.ma[1][2] = 1; A.ma[2][1] = 1; mat B=pow_mod(A, k); ans=(ans+B.ma[0][1]*max2+B.ma[0][2]*max1)%MOD; cout<<ans<<endl; } return 0; }
最后，只能参考学长的代码啦， 毕竟他山之石可以攻玉！。
#include<stdio.h> #include<string.h> #define mod 10000007 typedef long long ll; struct Mat { ll mat[5][5]; }; Mat operator * (Mat a, Mat b) { Mat c; memset(c.mat, 0, sizeof(c.mat)); for(int k=0; k<3; k++) for(int i=0; i<3; i++) for(int j=0; j<3; j++) c.mat[i][j] = (c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod; return c; } int main() { int n, k; Mat a, res; while(~scanf("%d%d", &n,&k)) { ll f1 = -1, f2 = -1, f, sum=0; for(int i=0; i<n; i++) { scanf("%lld", &f); sum += f; if(f > f1) f2 = f1, f1 = f; else if(f > f2) f2 = f; }//f1,f2为初始值 sum = sum-f1-f2; a.mat[0][0]=1; a.mat[0][1] = a.mat[0][2]=0; a.mat[1][0] = a.mat[1][1] = a.mat[1][2] = 1; a.mat[2][0] = a.mat[2][1] = 1; a.mat[2][2] = 0; for(int i=0; i<3; i++) for(int j=0; j<3; j++) res.mat[i][j] = (i==j); //初始为单位矩阵 for(; k>0; k>>=1) //快速幂 { if(k & 1) res = res*a; a = a*a; } sum += (f1+f2)*res.mat[0][0]+f1*res.mat[1][0]+f2*res.mat[2][0]; printf("%lld ", sum%mod); } return 0; }
感觉学长的代码和我的代码并没什么区别， 然而，，，，。 真是无语啦！。 不管怎样，弱渣仍需修练！
《6》这里还有几个特别巧妙地题。
http://www.cnblogs.com/kuangbin/category/405835.html（巨巨到底是巨巨！）