矩阵乘法和矩阵快速幂是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 }
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); }
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 }