矩陣乘法，矩陣快速冪學習筆記
矩陣乘法和矩陣快速冪是noip裡面經常考到的東西，可以用來優化一些dp，
例如菲波那契數列或者當前這一維於上一維所有狀態都有關係的dp
通常這些題目都有一個技巧，加入你的暴力dp式子推出來之後長得類似與f[ i ][ j ]=sum(f[ i-1 ][ p ]*g[ p ][ j ])
並且這個g陣列是一個固定的，並且大小不是很大的，可以預先處理出來的，那麼我們可以考慮運用矩陣快速冪進行優化。
那我們先從最基本的矩陣乘法開始入手，
我們定義3個矩陣 a，b，c;
a和b就是我們要進行想乘的兩個矩陣，c是我們的答案矩陣;
那麼這個矩陣乘法的狀態轉移式就是 c[ i ][ j ]=sum(a[ i ][ k ]*b[ k ][ j ]);
那麼貼個最普通也是用的範圍最廣的板子吧
1 struct matrix 2 { 3 int a[maxn][maxn]; 4 }x,y; 5 int n,m; 6 matrix cheng(matrix a,matrix b) 7 { 8 matrix c; 9 memset(c.a,0,sizeof(c.a)); 10 for(int i=1;i<=a矩陣第一維;i++) 11 { 12 for(int k=1;k<=a矩陣第二維和b矩陣第一維;k++) 13 { 14 for矩陣乘法(int j=1;j<=b矩陣的第二維;j++) 15 { 16 c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m; 17 } 18 } 19 } 20 return c; 21 }
那麼加入要乘的次數太多，肯定是要t的，所以我們可以利用一個普通的數字快速冪的做法，將數字相乘，換成矩陣相乘
程式碼實現其實跟普通的數字快速冪一樣的，只不過需要讓這個temp矩陣與a的第一次相乘結果為a，只需要把temp的整個第i維第i個賦值為1就好了
1 matrix qpow(matrix a,int n) 2 { 3 matrix b; 4 memset(b.a,0,sizeof(b.a)); 5 for(int i=1;i<=a矩陣最長的那一維;i++)//先定義一個矩陣，使他第一次乘等於他自己 6 { 7 b.a[i][i]=1; 8 } 9 while(n) 10 { 11 if(n&1) 12 { 13 b=cheng(b,a); 14 } 15 a=cheng(a,a); 16 n>>=1; 17 } 18 return b; 19 }快速冪
接下來講幾個題目吧：
1.菲波那契第N項
菲波那契的公式應該都知道吧就是f[ i ]=f[ i-1 ]+f[ i-2 ];
如果你一個一個來網上遞退，那麼肯定是會時空雙爆的，我們考慮矩陣快速冪進行優化;
但其實，我們是可以發現dp的規律的，定義f[ i ][ 1 ]為第f[ i ]個菲波那契，f[i][2]為i-1項，
那麼轉移式就是：f[ i ][ 1 ]=f[ i-1 ][ 1 ]+f[ i-1 ][ 2 ],f[ i ][ 2 ]=f[ i-1 ][1];
這貌似還沒有滿足矩陣快速冪的優化條件，因為還少一個數組，那麼我們應該如何做呢？
我們如果把加上變成加它乘以1，那麼這個表示式子就可以變成：
f[ i ][ 1 ]=f[ i-1 ][ 1 ]*1+f[ i-1 ][ 2 ]*1,f[ i ][ 2 ]=f[ i-1 ][1]×1+f[ i-1 ][ 2 ]*0;
我們可以根據這個定義一個數組：g[1][1]=1,g[1][2]=1,g[2][1]=1,g[2][2]=0;
再次轉移：f[ i ][ 1 ]=f[ i-1 ][ 1 ]*g[ 1 ][ 1 ]+f[ i-1 ][ 2 ]*g[ 2 ][ 1 ],
f[ i ][ 2 ]=f[ i-1 ][ 1 ]*g[ 1 ][ 2 ]+f[ i-1 ][ 2 ]*g[ 2 ][ 2 ];
發現規律了嗎？？？？
這不就是矩陣快速冪優化dp的通用公式嗎？？？？
那麼這個題就是f[1][1]=1，f[1][2]=1;的2×2的矩陣×g這個矩陣的（n-2）次冪，而最後的f[ 1 ][ 1 ]就是最後的答案！！
為啥是n-2次呢？
f[ 1 ][ 1 ]=1,f[ 1 ][ 2 ]=1時已經算到第2個了，你要算到第n個，那麼還需要算n-2次
程式碼？我都說到這分子上了，還需要我跟你貼程式碼？？？
還是貼一個吧
1 #include<bits/stdc++.h> 2 #define int long long 3 using namespace std; 4 inline int read() 5 { 6 int x=0,f=1; 7 char ch=getchar(); 8 while(ch<'0'||ch>'9') 9 { 10 if(ch=='-') 11 f=-1; 12 ch=getchar(); 13 } 14 while(ch>='0'&&ch<='9') 15 { 16 x=(x<<1)+(x<<3)+(ch^48); 17 ch=getchar(); 18 } 19 return x*f; 20 } 21 struct matrix 22 { 23 int a[3][3]; 24 }x,y; 25 int n,m; 26 matrix cheng(matrix a,matrix b) 27 { 28 matrix c; 29 memset(c.a,0,sizeof(c.a)); 30 for(int i=1;i<=2;i++) 31 { 32 for(int k=1;k<=2;k++) 33 { 34 for(int j=1;j<=2;j++) 35 { 36 c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m; 37 } 38 } 39 } 40 return c; 41 } 42 matrix qpow(matrix a,int n) 43 { 44 matrix b; 45 for(int i=1;i<=2;i++) 46 { 47 for(int j=1;j<=2;j++) 48 if(i==j) 49 b.a[i][j]=1; 50 else 51 b.a[i][j]=0; 52 } 53 while(n) 54 { 55 if(n&1) 56 { 57 b=cheng(b,a); 58 } 59 a=cheng(a,a); 60 n>>=1; 61 } 62 return b; 63 } 64 signed main() 65 { 66 n=read(),m=read(); 67 x.a[1][1]=1; 68 x.a[1][2]=1; 69 x.a[2][1]=0; 70 x.a[2][2]=0; 71 y.a[1][1]=y.a[2][1]=y.a[1][2]=1; 72 y.a[2][2]=0; 73 y=qpow(y,n-2); 74 x=cheng(x,y); 75 printf("%lld",x.a[1][1]); 76 }1
2.佳佳的菲波那契
給你一個新的定義s[ i ]=sum( f[ j ] )( 0<j<=i )
讓你求t[ i ]=sum( j*f[ j ] )( 0<j<=i )
菲波那契？這個我會，仔細看看？
這啥玩意，讓我求t還給我定義個s？
那麼這個題目肯定是有技巧的，當然，矩陣快速冪求菲波那契那都是基本操作了;
那麼我們考慮對s先進行化簡：s[n]=f[1]+----f[n]+1-1=f[1]+f[2]+f[2]----f[n]-1=f[n]+f[n+1]-1;
最後就是f[ n+2 ]-1,呀，這好像有點意思啊
那我們考慮對t進行化簡：
t[n]=n*f[n]+-------f[1];
t[n]=n*s[n]-(s[1]------s[n-1]);
t[n]=n*(f[n+2]-1)-(f[3]-----f[n+1]-(n-1))
t[n]=n*(f[n+2]-1)-(s[n+1]-(n-1)-2)
t[n]=n*f[n+2]-n-(f[n+3]-n-2)
t[n]=n*f[n+2]-f[n+3]+2
那麼最後直接出結果了好吧
得得得，看到這裡，如果你還不會做的話你可以扔掉你的腦子了
程式碼：
#include<bits/stdc++.h> #define int long long using namespace std; inline int read() { int x=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9') { if(ch=='-') f=-1; ch=getchar(); } while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); } return x*f; } struct matrix { int a[3][3]; }x,y; int n,m; matrix cheng(matrix a,matrix b) { matrix c; memset(c.a,0,sizeof(c.a)); for(int i=1;i<=2;i++) { for(int k=1;k<=2;k++) { for(int j=1;j<=2;j++) { c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m; } } } return c; } matrix qpow(matrix a,int n) { matrix b; for(int i=1;i<=2;i++) { for(int j=1;j<=2;j++) if(i==j) b.a[i][j]=1; else b.a[i][j]=0; } while(n) { if(n&1) { b=cheng(b,a); } a=cheng(a,a); n>>=1; } return b; } signed main() { n=read(),m=read(); x.a[1][1]=1; x.a[1][2]=1; x.a[2][1]=0; x.a[2][2]=0; y.a[1][1]=y.a[2][1]=y.a[1][2]=1; y.a[2][2]=0; y=qpow(y,n+1); x=cheng(x,y); int a=x.a[1][1],b=x.a[1][2]; printf("%lld",(b*n%m-a+2+m)%m); }2
3迷路scoi2009
這題，寧是否看出來如何使用矩陣快速冪？
如果寧翻到這裡了，那麼應該是看不出來了，那麼我來告訴你！！！！！
首先，我們先假設每個點的邊權是1那麼整個矩陣就是一個0 1矩陣對吧;
那麼我們可以想出一個dp暴力吧，
定義f[ i ][ j ][ k ]為花費k個體力從i到j的方案數，那麼狀態轉移方程就可以是
f[ i ][ j ][ k ]=sum(f[ i ][ p ][ k-1 ]*g[ p ][ j ])
我去，怎麼感覺有點眼熟啊。。。
已經很明顯了，我們把體力那一維度去掉，再仔細一考慮，這個矩陣k次方後的f[ i ][ j ]不就是i到j花費k個體力的方案數
哇，正解出來了！！！！！！！！！！
我要ac！！！！
10分！！！！！
不要著急，我們發現這個題目每一條邊不一定是1啊。。
那麼我們應該怎麼辦呢？？？？
假設i到j的體力為9，那麼是不是隻有k等9的時候，f[ i ][ j ]才能變成1;
那麼前8次怎麼辦呢？？？
我們發現，沒一條邊的邊權從0到9一共10中情況，那麼我們可以吧1個點分成10個點，組成一條邊權為1的鏈，
假設i到j是9，那麼我們可以讓i9與j1相連，
因為，i1與i2為1，i2與i3是1，所以這樣的話，k次方之後i1和j1就可以變成1了，
這道題目也就可以迎刃而解決了
最後f[ 1 ][ n*10-9 ]就是最後的答案！！！
程式碼
1 #include<bits/stdc++.h> 2 #define jb cout<<"jb"<<" "; 3 using namespace std; 4 inline int read() 5 { 6 int x=0,f=1; 7 char ch=getchar(); 8 while(ch<'0'||ch>'9') 9 { 10 if(ch=='-') 11 f=-1; 12 ch=getchar(); 13 } 14 while(ch>='0'&&ch<='9') 15 { 16 x=(x<<1)+(x<<3)+(ch^48); 17 ch=getchar(); 18 } 19 return x*f; 20 } 21 const int m=2009; 22 struct matrix 23 { 24 int a[201][201]; 25 void clear() 26 { 27 memset(a,0,sizeof(a));//不想再手寫memset了..... 28 } 29 }v; 30 int n,t,maxn; 31 inline matrix cheng(matrix a,matrix b) 32 { 33 matrix c; 34 c.clear(); 35 for(int i=1;i<=n;i++) 36 { 37 for(int k=1;k<=n;k++) 38 { 39 for(int j=1;j<=n;j++) 40 { 41 c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%m)%m; 42 } 43 } 44 } 45 return c; 46 } 47 inline matrix qpow(matrix a,int p) 48 { 49 matrix b; 50 b.clear(); 51 for(int i=1;i<=n;i++) 52 { 53 b.a[i][i]=1; 54 } 55 while(p) 56 { 57 if(p&1) 58 { 59 b=cheng(b,a); 60 } 61 a=cheng(a,a); 62 p>>=1; 63 } 64 return b; 65 } 66 inline int gd(int x,int y) 67 { 68 return (x-1)*10+y; 69 } 70 inline void init() 71 { 72 for(int i=1;i<=n;i++) 73 { 74 for(int j=1;j<=9;j++) 75 { 76 v.a[gd(i,j)][gd(i,j+1)]=1; 77 } 78 int c; 79 for(int j=1;j<=n;j++) 80 { 81 scanf("%1d",&c);//每次只吃掉1個數字，不寫字串型別的了 82 if(c) 83 { 84 v.a[gd(i,c)][gd(j,1)]=1; 85 } 86 } 87 } 88 } 89 int main() 90 { 91 n=read(),t=read(); 92 init(); 93 int N=n; 94 n*=10; 95 v=qpow(v,t); 96 printf("%d",v.a[1][10*N-9]); 97 }3