題目介紹：

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

#include<vector>
#include
 
<stdio.h>
#include<iostream>
#include<string.h>
using namespace std;
const int MAX_N = 220;
int E[MAX_N][MAX_N];
int M;
class Mat {
public:
    Mat(int N, int N2) {
        vc.resize(N);
        for (int i = 0; i < N; i++) {
            vc[i].resize(N2);
        }
    }
    vector
 
<int>& operator[](int idx) {
        return vc[idx];
    }
    int size() {
        return vc.size();
    }
private:
    vector<vector<int> > vc;
};
void multiply(Mat& CP, Mat& B, Mat& D) {
    memset(E, 0, sizeof(E));
    for (int i = 0; i < B.size(); i++) {
        for (int j = 0; j < B[0].size(); j++) {
            for (int k = 0; k < D.size(); k++) {
                // printf("prev, %d %d %d %d\n", E[i][j], B[i][k], D[k][j], B[i][k] * D[k][j]);
                E[i][j] = (E[i][j] + B[i][k] * (D[k][j] % M)) % M;
                // printf("last, %d %d %d %d\n", E[i][j], B[i][k], D[k][j], B[i][k] * D[k][j]);
            }
        }
    }
    for (int i = 0; i < B.size(); i++) {
        for (int j = 0; j < D[0].size(); j++) {
            CP[i][j] = E[i][j];
        }
    }
}

void debug(Mat& E, int N, int C) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < C; j++) {
            printf("%d ", E[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}
void mpow(Mat& B, int k) {
    int N = B.size();
    if (N == 0) {
        return;
    }
    int C = B[0].size();
    Mat cp(N, C);
    for (int i = 0; i < N; i++) {
        cp[i][i] = 1;
    }
    while (k) {
        if (k & 1) {
            // cp = cp * B^(2^k)
            multiply(cp, cp, B);
        }
        multiply(B, B, B);
        k = k >> 1;
    }
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            B[i][j] = cp[i][j];
        }
    }
}
int main() {
    int n, k;
    int A[MAX_N][MAX_N];
    while (cin >> n >> k >> M) {
        Mat B(n * 2, n * 2);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                cin >> A[i][j];
            }
        }
        // copy Mat to be:
        // | A | 0 |
        // | I | I |
        // 
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                B[i][j] = A[i][j];
            }
            B[n + i][i] = B[n + i][n + i] = 1;
        }
        mpow(B, k + 1);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                int a = B[n + i][j] % M;
                if (i == j) a = (a + M - 1) % M;
                printf("%d%c", a, j == n - 1 ? ‘\n‘ : ‘ ‘);
            }
        }
    }
    return 0;
}

【POJ 3233】矩陣乘積和 - 快速冪