POJ 3233 Matrix Power Series

矩阵类的题目。

这个k有点大，所以你单纯地快速幂然后求和是不可行的。

我们不妨令Sn=I+A+A^2+...+A^n

那么，Sn=Sn-1+An。写成矩阵的话我们有：

Sn =    I    I  *  Sn-1

An+1=0   A     An

此处的I是指单位矩阵。

然后矩阵快速幂即可。

#include <iostream>
#include <cstdio>
#include <cstring>

#define MAX 128

using namespace std;

typedef int i64;

i64 a[MAX][MAX], b[MAX][MAX], c[MAX][MAX], buff[MAX][MAX], vec[MAX], ans[MAX][MAX], res[MAX][MAX];

int MOD;

void matCpy(i64 a[MAX][MAX], i64 b[MAX][MAX], int n) {
    int i, j;
    
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            a[i][j] = b[i][j];
        }
    }
}

void norm(i64 a[MAX][MAX], int n) {
    int i, j;
    
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            a[i][j] = (i == j);
        }
    }
}

void matMul(const i64 a[MAX][MAX], const i64 b[MAX][MAX], i64 c[MAX][MAX], int n) {
    int i, j, k;
    
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            for (c[i][j] = k = 0; k < n; ++k) {
                c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % MOD;
            }
        }
    }
}

void matPow(i64 a[MAX][MAX], i64 b, i64 c[MAX][MAX], int n) {
    for (norm(c, n); b; b >>= 1) {
        if (b & 1) {
            matMul(c, a, buff, n);
            matCpy(c, buff, n);
        }
        matMul(a, a, buff, n);
        matCpy(a, buff, n);
    }
}

int main(){
    i64 n,k;
    while(scanf("%d%d%d",&n,&k,&MOD)!=EOF){
        if(n==0&&k==0)  break;
        i64 tem;
        memset(a,0,sizeof(a));
        memset(ans,0,sizeof(ans));
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                scanf("%d",&tem);
                a[n+i][n+j]=res[i][j]=tem%MOD;
            }
        }
        for(int i=0;i<n;i++){
                a[i][i]=a[i][n+i]=1;
        }
        /*
        for(int i=0;i<2*n;i++){
            for(int j=0;j<2*n;j++)  printf("%d ",a[i][j]);
            printf("\n");
        }*/
        matPow(a,k,b,2*n);
        /*
        for(int i=0;i<2*n;i++){
            for(int j=0;j<2*n;j++)  printf("%d ",b[i][j]);
            printf("\n");
        }
        */
        norm(c,n);
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                for(int k=0;k<n;k++){
                    ans[i][j]=(ans[i][j]+b[i][n+k]*res[k][j])%MOD;
                    ans[i][j]=(ans[i][j]+b[i][k]*c[k][j])%MOD;
                }
            }
        }
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                if(i==j)  printf("%d",(ans[i][j]-1+MOD)%MOD);
                else  printf("%d",ans[i][j]);
                if(j<n-1)  printf(" ");
            }
            printf("\n");
        }
        //printf("\n");
    }
    return 0;
}