poj 3233 Matrix Power Series

// 计算(A + A^2 +....+A^k) mod m ,其中 A为n*n矩阵

#include <iostream>        //利用转移矩阵和二分法
using namespace std; 
int n, m; 
struct node
{
    int matrix[31][62];    //把matrix拆分出两个n*n矩阵A,B,n最大为30,构造矩阵 E（A，B，0，1）,记录成E(A,B)
};    
node multiply(node a,node b)    //矩阵乘法 (A,B)* (C,D) =（A*C，A*D+B）
{
    node res;
    memset(res.matrix,0,sizeof(res.matrix));    //要显式初始化
    int i,j,k;
    for(i=0;i<n;++i)    //计算左上角矩阵A*C
        for(j=0;j<n;++j)
        {
            int& val=res.matrix[i][j];
            for(k=0;k<n;++k)        
                val=(val + a.matrix[i][k] * b.matrix[k][j])%m;    //要记得取余,由于m<10^4,所以两数相乘不会超int范围
        }
    for(i=0;i<n;++i)    //计算右上角矩阵A*D+B
        for(j=n;j<2*n;++j)
        {
            int& val=res.matrix[i][j];
            val=a.matrix[i][j];
            for(k=0;k<n;++k)        
                val=(val + a.matrix[i][k] * b.matrix[k][j])%m;
        }
    return res;
}
node pow(node p,int k)    //二分法计算矩阵p的k次方
{
    if(k==1)
        return p;
    else 
    {
        node q=pow(p,k/2);
        if(k%2==0)
            return multiply(q,q);
        else 
            return multiply( multiply(q,q) , p );
    }
}
int main()
{
    int k,i,j;
    node B;
    cin>>n>>k>>m;
    for(i=0;i<n;++i)    //左上角是原矩阵
        for(j=0;j<n;++j)
            cin>>B.matrix[i][j];
    for(i=0;i<n;++i)        //右上角是单位矩阵,主对角线为 1
        for(j=n;j<2*n;++j)
        {
            if(i==j-n)    //注意不是i==j
                B.matrix[i][j]=1;
            else
                B.matrix[i][j]=0;
        }

    B=pow(B,k+1);    //求B^(k+1)

    for(i=0;i<n;++i)        //B^(k+1)的右上角矩阵计算结果是1+A+A^2+.....+A^k,所以要减去 1 
        for(j=n;j<2*n;++j)
            if(i==j-n)
                B.matrix[i][j]--;
    for(i=0;i<n;++i)
    {
        for(j=n;j<2*n;++j)
        {
            cout<< ( B.matrix[i][j] + m ) % m <<" ";
        }
        cout<<endl;
    }
    return 0;
}


/*

给定A=n*n矩阵，计算A + A^2 +....+A^k的值。
                  A        1                              A^(k+1)      1+A+A^2+.....+A^k
构造矩阵 B   =                         ,易得B^(k+1) =                                                       

                  0         1                                 0                  1
(这里的1指的是单位矩阵，就是主对角线是1的单位矩阵)
所以问题转化为求B^(k+1)。
计算可知，在B(A,1,0,1)中，（0，1）不影响（A，1）的值，所以矩阵B可以简化为B（A，1）。
假设E(A,B)是E（A，B，0，1）的简化，F(C,D)是F（C，D，0，1）的简化，则E*F，矩阵的
乘法可定义为 E(A,B)*F(C,D) =（A*C，A*D+B，0，1），简化成 (A,B)* (C,D) =（A*C，A*D+B）
在求B^(k+1)时，可以利用二分法


*/