矩阵十题（3）

经典题目3
　　POJ3233
　　题目大意：给定矩阵A，求A + A^2 + A^3 + ... + A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据mod m。k<=10^9。
　　这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：
　　A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)
　　应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A + A^2 + A^3，即可得到原问题的答案。

奇数: F[n]=F[n-1]+A^n

偶数: F[n]=F[n/2]+F[n/2]*A^n/2

算法：二分+矩阵快速幂

代码如下：

 1 #include<stdio.h>
 2 #include<string.h>
 3 #define N 31
 4 struct Matrix
 5 {
 6     int a[N][N];
 7 }res,origin,tmp,A,B,C,ans;
 8 int n,m;
 9 Matrix mul(Matrix x,Matrix y)   //矩阵乘法
10 {
11     int i,j,k;
12     memset(tmp.a,0,sizeof(tmp.a));
13     for(i=1;i<=n;i++)
14         for(j=1;j<=n;j++)
15             for(k=1;k<=n;k++)
16             {
17                 tmp.a[i][j]+=(x.a[i][k]*y.a[k][j]);
18                 tmp.a[i][j]%=m;
19             }
20     return tmp;
21 }
22 Matrix quickpow(Matrix origin,int k)   //矩阵快速幂
23 {
24     int i;
25     memset(res.a,0,sizeof(res.a));
26     for(i=1;i<=n;i++)
27         res.a[i][i]=1;
28     while(k)
29     {
30         if(k&1)
31             res=mul(res,origin);
32         origin=mul(origin,origin);
33         k>>=1;
34     }
35     return res;
36 }
37 Matrix sum(Matrix x,Matrix y)   //矩阵求和
38 {
39     int i,j;
40     for(i=1;i<=n;i++)
41         for(j=1;j<=n;j++)
42             tmp.a[i][j]=(x.a[i][j]+y.a[i][j])%m;
43     return tmp;
44 }
45 Matrix bin(int k)
46 {
47     if(k<=1)
48         return A;
49     else if(k%2)   //奇数：F[n]=F[n-1]+A^n   
50     {
51         B=bin(k-1);
52         C=quickpow(A,k);   
53         return sum(B,C);
54     }
55     else     //偶数： F[n]=F[n/2]+F[n/2]*A^(n/2)
56     {
57         B=bin(k/2);
58         C=quickpow(A,k/2);
59         C=mul(C,B);
60         return sum(B,C);
61     }
62 }
63 int main()
64 {
65     int i,j,k;
66     scanf("%d%d%d",&n,&k,&m);
67     for(i=1;i<=n;i++)
68         for(j=1;j<=n;j++)
69             scanf("%d",&A.a[i][j]);
70     ans=bin(k);   //二分
71     for(i=1;i<=n;i++)
72     {
73         for(j=1;j<=n;j++)
74             if(j<n)
75                 printf("%d ",ans.a[i][j]);
76             else
77                 printf("%d\n",ans.a[i][j]);
78     }
79     return 0;
80 }

View Code