POJ 3233 Matrix Power Series（二分 / 矩阵套矩阵）

题目：http://poj.org/problem?id=3233

计算那个公式即可。

1.二分   偶数：

matrix b=solve(a,k/2);
return quick_pow(a,k/2)+1)+b;

奇数：

matrix b=solve(a,k/2);
return quick_pow(a,k)+quick_pow(a,(k-1)/2)*b+b;

2.矩阵套矩阵：

比起二分，我更喜欢这种思维的锻炼，但确实对我而言比较难写。

 A    1         A^2    1+A          A^3    1+A+A^2

 0   1            0         1             0             1

构造第一个矩阵，其他幂运算后去大矩阵的【0】【1】元素，减去单位矩阵

二分无限超时，而且我就是想用重载运算符，所以每次做很多多余运算，而第二种计算尽情用memset都能过（2844ms勉强过）

#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;

typedef int ll;
const int maxn=33;
int n,mod;
//定义矩阵乘法
struct matrix{
    ll arr[maxn][maxn];
    matrix operator*(matrix b){
        matrix ans;
        ll tmp;
        for(int i=0; i<maxn; i++)
        for(int j=0; j<maxn; j++){
            ans.arr[i][j] = 0;
            for(int k=0; k<maxn; k++){
                tmp = (arr[i][k]*b.arr[k][j])%mod;
                ans.arr[i][j] = (ans.arr[i][j] + tmp)%mod;
            }
        }
        return ans;
    }
    matrix operator+(matrix b){
        matrix ans;
        for(int i=0; i<maxn; i++)
        for(int j=0; j<maxn; j++){
            ans.arr[i][j]=(arr[i][j]+b.arr[i][j])%mod;
        }
        return ans;
    }
};
struct Matrix{
    matrix mat[2][2];
    Matrix operator*(Matrix b){
        Matrix ans;
        for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        for(int k=0;k<maxn;k++)
        for(int l=0;l<maxn;l++)
            ans.mat[i][j].arr[k][l]=0;
        for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        for(int k=0;k<2;k++){
            ans.mat[i][j]=ans.mat[i][j]+(mat[i][k]*b.mat[k][j]);
        }
        return ans;
    }
};

//矩阵快速幂
Matrix quick_pow(Matrix a,ll N){
    Matrix ans;
    for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
        memset(ans.mat[i][j].arr,0,sizeof(ans.mat[i][j].arr)); //开到了maxn，所以后面也要初始化
    for(int i=0; i<n; i++)
        ans.mat[0][0].arr[i][i]=ans.mat[1][1].arr[i][i]=1;
    while(N){
        if(N&1)
            ans = ans*a;
        a = a*a;
        N >>= 1;
    }
    return ans;
}
int main(){
    int k;
    Matrix MM;
    scanf("%d%d%d",&n,&k,&mod);
    for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
        memset(MM.mat[i][j].arr,0,sizeof(MM.mat[i][j].arr));
    for(int i=0;i<n;i++){
        MM.mat[0][1].arr[i][i]=MM.mat[1][1].arr[i][i]=1;
        for(int j=0;j<n;j++){
            scanf("%d",&MM.mat[0][0].arr[i][j]);
            MM.mat[1][0].arr[i][j]=0;
            if(i!=j)
                MM.mat[0][1].arr[i][j]=MM.mat[1][1].arr[i][j]=0;
        }
    }
    Matrix ans=quick_pow(MM,k+1);
    for(int i=0;i<n;i++){
        ans.mat[0][1].arr[i][i]--;
        if(ans.mat[0][1].arr[i][i]<0)
            ans.mat[0][1].arr[i][i]+=mod;
        printf("%d",ans.mat[0][1].arr[i][0]);
        for(int j=1;j<n;j++){
            printf(" %d",ans.mat[0][1].arr[i][j]);
        }
        printf("\n");
    }
    return 0;
}