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 + 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

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[N][N];
	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;
}