矩陣乘法是一種高效的算法可以把一些一維遞推優化到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\n", 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\n", 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\n", 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\n",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)?'\n':' '); } } } 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\n", 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\n", 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\n", sum%mod); } return 0; }
感覺學長的代碼和我的代碼并沒什么區別， 然而，，，，。 真是無語啦！。 不管怎樣，弱渣仍需修練！
《6》這里還有幾個特別巧妙地題。
http://www.cnblogs.com/kuangbin/category/405835.html（巨巨到底是巨巨！）
文章列表
留言列表
