POJ 3233 Matrix Power Series 【矩阵快速幂，矩阵加速】

| A+A^2+A^3+…Ak |   |A   A| (k-1）次方   | A |
|                | = |     |              |   |
|     E          |   |0   E|              | E |

A是输入的矩阵

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

#define N 66
int n, a[N][N], Mod;

//c = a*b
void Multi(int a[][N], int b[][N], int c[][N]) {
    for (int i=0; i<n; i++)
        for (int j=0; j<n; j++) {
            c[i][j] = 0;
            for (int k=0; k<n; k++)
                c[i][j] = (c[i][j] + a[i][k]*b[k][j]) % Mod;
        }
}
//d = s
void copy(int d[][N], int s[][N]) {
    for (int i=0; i<n; i++) for (int j=0; j<n; j++)
        d[i][j] = s[i][j];
}
//a = a^k % Mod
void PowerMod(int a[][N], int b) {
    int t[N][N], ret[N][N];
    for (int i=0; i<n; i++) ret[i][i] = 1;
    while (b) {
        if (b & 1) { Multi(ret, a, t); copy(ret, t); }
        Multi(a, a, t); copy(a, t);
        b >>= 1;
    }
    copy(a, ret);
}
void print(int x[][N], int y) {
    for (int i=0; i<y; i++) {
        printf("%d", x[i][0]);
        for (int j=1; j<y; j++) printf(" %d", x[i][j]);
        printf("\n");
    }
}
int main() {
    int k;
    scanf("%d%d%d", &n, &k, &Mod);
    memset(a, 0, sizeof(a));

    for (int i=0; i<n; i++) for (int j=0; j<n; j++) scanf("%d", &a[i][j]);
    if (k > 1) {
        int m = n, b[N][N], c[N][N];
        memset(b, 0, sizeof(b));
        for (int i=0; i<n; i++) for (int j=0; j<n; j++) {
            a[i][j+n] = a[i][j];
            a[i+n][i+n] = 1;
            b[i][j] = a[i][j];
            b[i+n][i] = 1;
        }
        n = n*2;
        PowerMod(a, k-1);
        memset(c, 0, sizeof(c));
        for (int i=0; i<n; i++)
            for (int j=0; j<m; j++)
                for (int k=0; k<n; k++)
                    c[i][j] = (c[i][j] + a[i][k]*b[k][j]) % Mod;
        print(c, n/2);
    } else print(a, n);
    return 0;
}