POJ3233 Matrix Power Series (等比矩阵求和)
题意::求等比矩阵和:S = A + A2 + A3 + … + Ak
参考巨佬们的经典算法( 接近模板 )
:::构造矩阵 E为单位矩阵
看不明白的可以手模一下,
答案就是 右上角子矩阵减一个同阶单位矩阵
1 #include<bits/stdc++.h> 2 #define ll long long 3 using namespace std; 4 const int maxn=65; 5 6 int n,MOD,k; 7 struct mat 8 { 9 int a[maxn][maxn]; 10 mat(){ 11 memset(a,0,sizeof(a)); 12 } 13 }; 14 mat mul(mat x,mat y) 15 { 16 mat res; 17 for(int i=0;i<n;i++) 18 { 19 for(int j=0;j<n;j++) 20 { 21 for(int k=0;k<n;k++) 22 { 23 res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j]%MOD)%MOD; 24 } 25 } 26 } 27 return res; 28 } 29 mat qpow(mat x,int p) 30 { 31 mat res; 32 for(int i=0;i<n;i++){ 33 res.a[i][i]=1; 34 } 35 while(p) 36 { 37 if(p&1){ 38 res=mul(res,x); 39 } 40 x=mul(x,x); 41 p>>=1; 42 } 43 return res; 44 } 45 int main() 46 { 47 scanf("%d%d%d",&n,&k,&MOD); 48 mat ans; 49 for(int i=0;i<n;i++){ 50 for(int j=0;j<n;j++){ 51 scanf("%d",&ans.a[i][j]); 52 ans.a[i][j]%=MOD; 53 } 54 } 55 for(int i=n;i<(n<<1);i++){ 56 ans.a[i-n][i]=ans.a[i][i]=1; 57 } 58 n=n<<1; 59 ans=qpow(ans,k+1); 60 n=n>>1; 61 for(int i=0;i<n;i++) 62 { 63 ans.a[i][i+n]=(ans.a[i][i+n]-1+MOD)%MOD; 64 } 65 for(int i=0;i<n;i++) 66 { 67 for(int j=n;j<(n<<1);j++) 68 { 69 printf("%d ",ans.a[i][j]); 70 } 71 printf("\n"); 72 } 73 return 0; 74 }
矩阵快速幂+二分
S = A + A2 + A3 + … + Ak
对于整数而言::
n为偶数时::S[n]=(a^(n/2)+1)*S[n/2]
n为奇数时::S[n]=(a^(n/2+1)+1)*S[n/2]+(a^(n/2)+1);
对于矩阵也同样适用 ,当n==0时 S为单位矩阵;
1 #include<bits/stdc++.h> 2 #define ll long long 3 using namespace std; 4 const int maxn=65; 5 6 int n,MOD,k; 7 struct mat 8 { 9 int a[maxn][maxn]; 10 mat(){ 11 memset(a,0,sizeof(a)); 12 } 13 }; 14 mat mul(mat x,mat y) 15 { 16 mat res; 17 for(int i=0;i<n;i++){ 18 for(int j=0;j<n;j++){ 19 for(int k=0;k<n;k++){ 20 res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j]%MOD)%MOD; 21 } 22 } 23 }return res; 24 } 25 mat add(mat x,mat y) 26 { 27 mat res; 28 for(int i=0;i<n;i++){ 29 for(int j=0;j<n;j++){ 30 res.a[i][j]=(res.a[i][j]+x.a[i][j]+y.a[i][j])%MOD; 31 } 32 } 33 return res; 34 } 35 mat qpow(mat x,int p) 36 { 37 mat res; 38 for(int i=0;i<n;i++){ 39 res.a[i][i]=1; 40 } 41 while(p){ 42 if(p&1){ 43 res=mul(res,x); 44 } 45 x=mul(x,x); 46 p>>=1; 47 } 48 return res; 49 } 50 mat solve(mat x,int k) 51 { 52 if(k==1){ 53 return x; 54 } 55 mat t=solve(x,k/2); 56 if(k&1){ 57 mat res=qpow(x,k/2+1); 58 t=add(add(mul(res,t),t),res); 59 } 60 else{ 61 mat res=qpow(x,k/2); 62 t=add(mul(res,t),t); 63 } 64 return t; 65 } 66 67 int main() 68 { 69 scanf("%d%d%d",&n,&k,&MOD); 70 mat ans; 71 for(int i=0;i<n;i++){ 72 for(int j=0;j<n;j++){ 73 scanf("%d",&ans.a[i][j]); 74 ans.a[i][j]%=MOD; 75 } 76 } 77 ans=solve(ans,k); 78 for(int i=0;i<n;i++){ 79 for(int j=0;j<n;j++){ 80 printf("%d ",ans.a[i][j]); 81 }printf("\n"); 82 } 83 return 0; 84 }
纵使单枪匹马,也要勇闯天涯