POJ 3233 Matrix Power Series

　　题目的意思是已知一个n阶矩阵A，以及K，M，求S = A + A^2 + A^3 + ... + A^K，由于里面的数值当K很大时偏大，对M取模。其中数据的范围为：

　  n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 

　　构造分块矩阵 B = A  I

　　　　　　　　　　  0  I

　　由数学归纳法可证得 B^(K+1）= A^K  I+A^1+A^2+...+A^K

　　　　　　　　　　　　　　　　　　0     I

　　这样这道题目就纯粹转化为了求B矩阵的K+1次方了，二分即可。

　　代码如下：

#include <iostream>

#define MAXN 61

using namespace std;

int NN , N , M , K , temp;
int a[MAXN][MAXN] , ans[MAXN][MAXN];
int b[32][MAXN][MAXN] , u[31];

void    init(int T) {
    memset(b,0,sizeof(b));
    for (int i = 0;i < NN;i++)
        for (int j = 0;j < NN;j++)
            b[1][i][j] = a[i][j];
    u[1] = 1; temp = 1;
    while (u[temp]*2 < T) {
        temp ++;
        u[temp] = u[temp-1] * 2;
        for (int k = 0;k < NN;k++)
            for (int i = 0;i < NN;i++)
                for (int j = 0;j < NN;j++) {
                    b[temp][i][j] += (b[temp-1][i][k] * b[temp-1][k][j]) % M;
                    b[temp][i][j] %= M;
                }
    }
    memset(ans,0,sizeof(ans));
    for (int i = 0;i < NN;i++)  ans[i][i] = 1;
}

void    calc(int T) {
    int bak[MAXN][MAXN];
    while (T) {
        while (T < u[temp]) temp --;
        memcpy(bak,ans,sizeof(bak));
        memset(ans,0,sizeof(ans));
        for (int k = 0;k < NN;k++)
            for (int i = 0;i < NN;i++)
                for (int j = 0;j < NN;j++) {
                    ans[i][j] += (bak[i][k] * b[temp][k][j]) % M;
                    ans[i][j] %= M;
                }
        T -= u[temp];
    }
    for (int i = 0;i < NN;i++)
      ans[i][i+N] = (ans[i][i+N] - 1 + M) % M;
}

int main () {
    memset(a,0,sizeof(a));
    cin >> N >> K >> M;
    for (int i = 0;i < N;i++)
        for (int j = 0;j < N;j++)
            cin >> a[i][j];
    NN = N*2;
    for (int i = 0;i < N;i++)
        a[i][i+N] = 1;
    for (int i = N;i < NN;i++)
        a[i][i] = 1;
    memcpy(ans,a,sizeof(a));
    init(K+1);
    calc(K+1);
    for (int i = 0;i < N;i++) {
        for (int j = N;j < NN-1;j++)
            cout << ans[i][j] << " ";
        cout << ans[i][NN-1] << endl;
    }
    return 0;
}

　　测试结果如下：