poj 3233Matrix Power Series（矩阵快速幂 二分求和 求累乘的和）

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 20868		Accepted: 8732

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). 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

Source

POJ Monthly--2007.06.03, Huang, Jinsong

/*
题目大意：给你A矩阵，A矩阵是n*n的一个矩阵，现在要你求S = A + A^2 + A^3 + … + A^k.
那么s一定也是一个N*N的矩阵，最后要你输出s，并且s的每一个元素对m取余数

解题思路：因为S可以看成S=A(I+A(I+A(I+...A(I+A)))) (I是单位矩阵)
拿k=3举例S=A(I+A(I+A))
那么我们想，可不可以构造一个矩阵T使得T*T（因为是k次幂）这样乘下去每次可以得到A*(A+I)
那么肯定T有个两个元素就是A与I
那么假设：T={A  I }
I  I
那么T=T*T={A*A+I*I       A*I+I*I}
A*I+I*I        I*I+I*I
这样存在一个I*（A+I）的式子 ，当T再乘以T的时候会出现A（A+I）
这个时候我们可以简化将T={A  I}
0   I
这样可以简化很多计算T*T={A*A   A*I+I*I}
0           I
那么容易得到T^(K+1)={A^(K+1)           I+A+A^2+A^3+...+A^K}
0                                I
这样我们只需要算T的k+1次幂就可以了

而k如此庞大所以需要二分来对T求k+1次幂
对T求快速幂：
首先我们知道A^19=A^16*A^2*A^1，因为19=B(10011)

那么就这样，拿A总是去和本身去乘，那么就可以取到 A A^2 A^4 A^8 A^16...

然后问A A^2 A^4 A^8 A^16 ...这么多项取与不取的抉择
其实就是一个求二进制的这样一个过程
19%2=1 那么A取进来
19/2=9
9%2=1 那么A^2取进来
9/2=4
4%2=0 那么A^4不用取进来
4/2=2
2%2=0 那么A^8不用去进来
2/2=1
1%2=1 那么A^16取进来
1/2=0 计算完毕
*/

#include<iostream>
#include<string>
#include<cstring>
#include<cstdio>
using namespace std;
#define N 70
struct matrix{
int mat

;
matrix()
{
memset(mat,0,sizeof(mat));
}

};
matrix d;
int n,m,k;
int Mod;
matrix mul(matrix a,matrix b)
{
matrix c;
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
c.mat[i][j]=0;
for(int k=0;k<n;k++)
{
c.mat[i][j]+=(a.mat[i][k]*b.mat[k][j]);
if(c.mat[i][j]>=Mod)
c.mat[i][j]%=Mod;
}
c.mat[i][j]%=Mod;
}
}
return c;
}
matrix expo(matrix a,int k)
{
if(k==1)
return a;
matrix e;
for(int i=0;i<n;i++)
{
e.mat[i][i]=1;
}
if(k==0)
return e;
while(k)
{
if(k&1)
e=mul(a,e);
a=mul(a,a);
k>>=1;
}
return e;
}
matrix add(matrix a,matrix b)
{
matrix t;
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
t.mat[i][j]=(a.mat[i][j]+b.mat[i][j]);
if(t.mat[i][j]>=Mod)
t.mat[i][j]%=Mod;
}
}
return t;
}
matrix sum(int k)//类似于树状数组求和  sum(4)包含了e e^2 e^3 e^4
{
if(k==1)
{
return d;
}
if(k&1)
{
return add(sum(k-1),expo(d,k));    //以19为例产生递归   sum18+e^19  sum9+e^9  sum8+e^9  sum4*e^4  sum2*e^2  sum1*e  e  最终就是累加到19
}
else {
matrix s=sum(k>>1);
return add(s,mul(s,expo(d,k>>1)));
}
}
int main()
{
while(~scanf("%d%d%d",&n,&k,&m))
{
matrix ans,t;
Mod=m;
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
scanf("%d",&d.mat[i][j]);
if(d.mat[i][j]>=m)
{
d.mat[i][j]%=m;
}
}
}
ans=sum(k);
print(ans);
}
return 0;
}