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;
}