POJ3233 Matrix Power Series (等比矩阵求和)

题意：：求等比矩阵和：S = A + A² + A³ + … + A^k

参考巨佬们的经典算法（接近模板）

:::构造矩阵 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("
");
72     }
73     return 0;
74 }

View Code

矩阵快速幂+二分

S = A + A² + A³ + … + A^k

对于整数而言：：

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("
");
82     }
83     return 0;
84 }

View Code

纵使单枪匹马，也要勇闯天涯