POJ 3233 矩阵快速幂 等比数列

//http://www.cnblogs.com/jiangjing/archive/2013/05/28/3103336.html

/*
分析：求a^1+..a^n这是矩阵乘法中关于等比矩阵的求法：

|A  E|

|0  E|

其中的A为m阶矩阵，E是单位矩阵，0是零矩阵。而我们要求的是：

A^1+A^2+..A^L，由等比矩阵的性质

|A  ,  1|                 |A^n , 1+A^1+A^2+....+A^(n-1)|

|0  ,  1| 的n次方等于     |0   ,         1             |

所以我们只需要将A矩阵扩大四倍，变成如上形式的矩阵B，然后开L+1次方就可以得到1+A^1+A^2+....+A^L。由于多了一个1，所以最后得到的答案我们还要减去1。同理我们把矩阵A变成B：

|A  E|

|0  E|

然后我们就是求B的n+1次幂之后得到的矩阵为|x1   x2|

|x3   x4|

右上角的矩阵x2减去单位矩阵E，得到就是要求的矩阵了！

*/

#include<stdio.h>
#include<string>
#include<cstring>
#include<queue>
#include<algorithm>
#include<functional>
#include<vector>
#include<iomanip>
#include<math.h>
#include<iostream>
#include<sstream>
#include<stack>
#include<set>
#include<bitset>
using namespace std;
const int INF=0x3f3f3f3f;
const int MOD=10000;
const int SIZE=64;
typedef long long ll;
struct Matrix
{
int n;
int Mod;
ll Mat[SIZE][SIZE];
Matrix(int a):n(a),Mod(MOD)
{
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
Mat[i][j]=0;
}
Matrix(int a,int m):n(a),Mod(m)
{
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
Mat[i][j]=0;
}
Matrix operator + (Matrix x)
{
Matrix result(n,Mod);
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result.Mat[i][j]=(Mat[i][j]%Mod+x.Mat[i][j]%Mod+Mod)%Mod;
return result;
}
Matrix operator * (Matrix x)
{
Matrix result(n,Mod);
for (int k=0; k<n; k++)
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result.Mat[i][j]=(result.Mat[i][j]+Mat[i][k]%Mod*x.Mat[k][j]%Mod+Mod)%Mod;
return result;
}
Matrix operator ^ (ll x)
{
Matrix temp(n,Mod),a(*this);
for (int i=0; i<n; i++)
temp.Mat[i][i]=1;
while (x)
{
if (x&1)
temp=a*temp;
a=a*a;
x>>=1;
}
return temp;
}
void print()
{
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
cout<<Mat[i][j]<<(j==n-1?'\n':' ');
}
};
int main()
{
cin.sync_with_stdio(false);
int N,K,M;
cin>>N>>K>>M;
Matrix Ans(2*N,M);
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
cin>>Ans.Mat[i][j];
for (int i=0; i<N; i++)
Ans.Mat[i][i+N]=1;
for (int i=N; i<2*N; i++)
Ans.Mat[i][i]=1;
Ans=Ans^(K+1);
for (int i=0; i<N; i++)
Ans.Mat[i][i+N]=(Ans.Mat[i][i+N]-1+Ans.Mod)%Ans.Mod;
for (int i=0; i<N; i++)
for (int j=N; j<2*N; j++)
cout<<Ans.Mat[i][j]<<(j==2*N-1?'\n':' ');
return 0;
}