POJ 3233 Matrix Power Series (矩阵快速幂+二分求解)
题意:求S=(A+A^2+A^3+...+A^k)%m的和
方法一:二分求解
S=A+A^2+...+A^k
若k为奇数:
S=(A+A^2+...+A^(k/2))+A^(k/2)*(A+A^2+...+A^(k/2))+A^k
若k为偶数:
S=(A+A^2+...+A^(k/2))+A^(k/2)*(A+A^2+...+A^(k/2))
也可以这么二分(其实和前面的差不多):
S(2n)=A+A^2+...+A^2n=(1+A^n)*(A+A^2+...+A^n)=(1+A^n)*S(n)
S(2n+1)=A+A^2+...+A^(2n+1)=A(1+A+..+A^2n)=A+(A+A^(n+1))*S(n)
一开始1900+ms,优化了下1500ms...还是太慢了。。。
本来在递归的时候,用快速幂计算A^(k/2)
后来改用直接递归的同时,计算A^(k/2),直接变成200ms左右。。。瞬间提升了10倍。。。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include <iostream> #include <cstdio> #include <string.h> using namespace std; const int maxn=31; int mod; int n,k,m; struct Matrix{ int m[maxn][maxn]; void init(){ memset(m,0,sizeof(m)); } void initE(){ memset(m,0,sizeof(m)); for(int i=0;i<n;i++) m[i][i]=1; } }A; //重载+运算符 Matrix operator+(Matrix a,Matrix b){ Matrix c; for(int i=0;i<n;i++){ for(int j=0;j<n;j++) c.m[i][j]=(a.m[i][j]+b.m[i][j])%mod; } return c; } //重载*运算符 Matrix operator*(Matrix a,Matrix b){ Matrix c; for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ c.m[i][j]=0; for(int k=0;k<n;k++){ c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod; } } } return c; } //矩阵快速幂 Matrix MquickPow(Matrix A,int b){ Matrix ret; ret.initE(); while(b){ if(b&1) ret=ret*A; A=A*A; b=b>>1; } return ret; } Matrix p; Matrix dfs(Matrix A,int k){ if(k==1){ p=A; return A; } Matrix ret,ans; ret=dfs(A,k/2); //Matrix p=MquickPow(A,k/2);如果用快速幂计算p=A^(k/2),则要1500ms,而直接在递归的时候同时计算p,则只要188ms。 if(k&1){ //return ret+ret*p+p*p*A; ans=ret+ret*p+p*p*A; p=p*p*A; } else{ //return ret+ret*p; ans=ret+ret*p; p=p*p; } return ans; } int main() { scanf("%d%d%d",&n,&k,&m); mod=m; for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ scanf("%d",&A.m[i][j]); } } Matrix ans; ans=dfs(A,k); for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ printf("%d ",ans.m[i][j]); } printf("\n"); } return 0; }
方法二:
688ms
http://blog.sina.com.cn/s/blog_9d5278a301015mbd.html
http://blog.csdn.net/wangjian8006/article/details/7868864
S=A(E+A(E+A(E+...A(E+A))))
这样就可以想,不妨构造一个矩阵T使得T*T,这样乘下去每次可以得到A*(A+E)+E
A 0 A^2 0 A^3 0
B= B^2= B^3=
E E A+E E A^2+A+E E
不难得出: A^(k+1) 0
B^(k+1)=
S(k) E
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include <iostream> #include <iostream> #include <cstdio> #include <string.h> using namespace std; const int maxn=65; int mod; int n,k,m; struct Matrix{ int m[maxn][maxn]; void init(){ memset(m,0,sizeof(m)); } void initE(){ memset(m,0,sizeof(m)); for(int i=0;i<n;i++) m[i][i]=1; } }A; //重载+运算符 Matrix operator+(Matrix a,Matrix b){ Matrix c; for(int i=0;i<n;i++){ for(int j=0;j<n;j++) c.m[i][j]=(a.m[i][j]+b.m[i][j])%mod; } return c; } //重载*运算符 Matrix operator*(Matrix a,Matrix b){ Matrix c; for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ c.m[i][j]=0; for(int k=0;k<n;k++){ c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod; } } } return c; } //矩阵快速幂 Matrix MquickPow(Matrix A,int b){ Matrix ret; ret.initE(); while(b){ if(b&1) ret=ret*A; A=A*A; b=b>>1; } return ret; } int main() { Matrix B; B.init(); scanf("%d%d%d",&n,&k,&m); mod=m; for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ scanf("%d",&B.m[i][j]); } B.m[i+n][i]=1; B.m[i+n][i+n]=1; } n=n*2; Matrix ans=MquickPow(B,k+1); for(int i=n/2;i<n;i++){ for(int j=0;j<n/2;j++){ if(i==j+n/2) ans.m[i][j]=(ans.m[i][j]-1+mod)%mod; //对角线还要减去单位矩阵的1 printf("%d ",ans.m[i][j]); } printf("\n"); } return 0; }