POJ 3233 Matrix Power Series(矩阵快速幂+二分求和)

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

题意：给定一个n*n的矩阵，求矩阵k阶阶乘，即求A^1+A^2+......+A^k的和
思路：因为k的范围是1到10^9，所以不能暴力去做，于是应该想到比较快的时间复杂度：二分。在快速幂的基础上，通过二分来求和

具体代码如下：

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cstring>
  4 #include <cstdlib>
  5 #include <cmath>
  6 #include <algorithm>
  7 using namespace std;
  8 
  9 const int m=33;
 10 
 11 struct node
 12 {
 13     int matrix[m][m];
 14     node(){}
 15     friend node operator * (node a,node b);
 16     friend node operator + (node a,node b);
 17 };
 18 
 19 int mod;
 20 int n;
 21 
 22 void init(node &a)
 23 {
 24     for(int i=0;i<n;i++)
 25     {
 26         for(int j=0;j<n;j++)
 27             a.matrix[i][j]=0;
 28     }
 29 
 30 }
 31 
 32 node operator * (node a,node b)
 33 {
 34     node c;
 35     init(c);
 36     for(int k=0;k<n;k++)
 37     {
 38         for(int i=0;i<n;i++)
 39         {
 40             for(int j=0;j<n;j++)
 41             {
 42                 c.matrix[i][j]=(c.matrix[i][j]+a.matrix[i][k]*b.matrix[k][j])%mod;
 43             }
 44         }
 45     }
 46     return c;
 47 }
 48 
 49 node operator + (node a,node b)
 50 {
 51     node tmp;
 52     init(tmp);
 53     for(int i=0;i<n;i++)
 54     {
 55         for(int j=0;j<n;j++)
 56             tmp.matrix[i][j]=(a.matrix[i][j]+b.matrix[i][j])%mod;
 57     }
 58     return tmp;
 59 }
 60 
 61 node fast_pow(node a,long long k)
 62 {
 63     node b;
 64     for(int i=0;i<n;i++)
 65     {
 66         for(int j=0;j<n;j++)
 67         {
 68             if(i==j) b.matrix[i][j]=1;
 69             else b.matrix[i][j]=0;
 70         }
 71     }
 72     while(k)
 73     {
 74         if(k&1) b=a*b;
 75         a=a*a;
 76         k>>=1;
 77     }
 78     return b;
 79 }
 80 
 81 node sum(node a,long long k)    //二分求和
 82 {
 83     node tmp1,tmp2;
 84     if(k==1) return a;
 85     tmp1=sum(a,k/2);      //递归
 86     if(k&1)
 87     {
 88         tmp2=fast_pow(a,k/2+1);
 89         tmp1=tmp2*tmp1+tmp1;
 90         return tmp1+tmp2;
 91     }
 92     else
 93     {
 94         tmp2=fast_pow(a,k/2);
 95         return tmp2*tmp1+tmp1;
 96     }
 97 }
 98 
 99 int main()
100 {
101     long long k;
102     while(~scanf("%d%I64d%d",&n,&k,&mod))
103     {
104         node a;
105         for(int i=0;i<n;i++)
106         {
107             for(int j=0;j<n;j++)
108             {
109                 scanf("%d",&a.matrix[i][j]);
110                 a.matrix[i][j]%= mod;
111             }
112         }
113         node b=sum(a,k);
114         for(int i=0;i<n;i++)
115         {
116             for(int j=0;j<n;j++)
117             {
118                 printf("%d%c",b.matrix[i][j],j==n-1?'
':' ');
119             }
120         }
121     }
122     return 0;
123 }