学步园

矩阵乘法。

S = A + A²
+ A³
+ … + A^k
，把问题转化以加速，令

B = A  I

      0  I

则B^(k + 1) = A^(k + 1)      I + A + A²
+ A³
+ … + A^k

                            0                          I

用二分法求B^(k + 1)

来源：http://archive.cnblogs.com/a/1960189/

#include<iostream><br /> #include<cstdio><br /> #include<cstring><br /> using namespace std;<br /> class node<br /> {<br /> public:<br /> int matrix[66][66];<br /> }ans,b;<br /> int n,m,k;<br /> void init(int n,int m)<br /> {<br /> int i,j;<br /> memset(ans.matrix ,0,sizeof(ans.matrix ));<br /> memset(b.matrix ,0,sizeof(b.matrix ));<br /> for(i=1;i<=n;i++)<br /> for(j=1;j<=n;j++)<br /> {<br /> scanf("%d",&ans.matrix [i][j]);<br /> ans.matrix [i][j]=ans.matrix [i][j]%m;<br /> b.matrix [i][j]=ans.matrix [i][j];<br /> }<br /> for(i=1;i<=n;i++) {ans.matrix [i][n+i]=1; b.matrix [i][n+i]=1;}<br /> for(i=1+n;i<=n*2;i++) { ans.matrix [i][i]=1;b.matrix [i][i]=1;}<br /> }<br /> void mul(node &a,node &b)<br /> {<br /> node x;<br /> int i,k,j,sum;<br /> for(i=1;i<=n*2;i++)<br /> for(j=1;j<=n*2;j++)<br /> x.matrix[i][j]=a.matrix[i][j];<br /> for(i=1;i<=n*2;i++)<br /> for(j=1;j<=n*2;j++)<br /> {<br /> sum=0;<br /> for(k=1;k<=n*2;k++)<br /> sum+=(x.matrix [i][k]*b.matrix [k][j]);<br /> a.matrix [i][j]=sum%m;<br /> }<br /> }<br /> void solve()<br /> {<br /> node a;<br /> int i,j;<br /> for(i=1;i<=n*2;i++)<br /> for(j=1;j<=n*2;j++)<br /> {<br /> a.matrix [i][j]=ans.matrix [i][j];<br /> }<br /> mul(ans,a);<br /> }<br /> void work(int num)<br /> {<br /> if(num==1) return ;<br /> work(num/2);<br /> solve();<br /> if(num%2==1)<br /> mul(ans,b);<br /> }<br /> int main()<br /> {<br /> int i,j;<br /> cin>>n>>k>>m;<br /> init(n,m);<br /> work(k+1);<br /> for(i=1;i<=n;i++)<br /> {<br /> if(ans.matrix [i][i+n]==0)<br /> ans.matrix [i][i+n]=m;<br /> ans.matrix [i][i+n]-=1;<br /> }<br /> for(i=1;i<=n;i++)<br /> {<br /> cout<<ans.matrix [i][n+1];<br /> for(j=2+n;j<=n*2;j++)<br /> cout<<" "<<ans.matrix [i][j];<br /> cout<<endl;<br /> }<br /> return 0;<br /> }