poj 3233(矩阵快速幂)

题意：给出一个矩阵A和数字k，要求出矩阵S = A + A^2 + A^3 + … + A^k。
题解：首先A^x可以计算，然后需要折半计算，比如s(k) = (1 + A^(k/2)) * s(k/2)，但k的奇偶不同需要分情况。

#include <stdio.h>
#include <string.h>
const int N = 35;
struct Mat {
    int g[N][N];
};
int n, k, m;

Mat mul(Mat x, Mat y) {
    Mat temp;
    memset(temp.g, 0, sizeof(temp.g));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                temp.g[i][j] = (temp.g[i][j] + x.g[i][k] * y.g[k][j]) % m;
    return temp;
}

Mat add(Mat x, Mat y) {
    Mat temp;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            temp.g[i][j] = (x.g[i][j] + y.g[i][j]) % m;
    return temp;
}

Mat Pow(Mat x, int k) { //计算 A^k
    if (k == 0) { //单位阵
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                x.g[i][j] = (i == j) ? 1 : 0;
        return x;
    }
    if (k == 1)
        return x;
    Mat temp = Pow(x, k / 2); // temp = A^(k/2)
    if (k % 2)
        return mul(mul(temp, temp), x);//exp: A^7 = A^3 * A^3 * A;
    return mul(temp, temp);
}

Mat Cal(Mat x, int k) { //计算 s(k) = A + A^2 + A^3 + ... + A^k
    if (k == 1)
        return x;
    Mat a = Pow(x, (k + 1) / 2); // a = A^((k+1)/2)
    Mat b = Cal(x, k / 2); // b = s(k/2)
    if (k % 2 == 0)
        return mul(add(Pow(x, 0), a), b); //exp: s(6) = (1 + A^3) * s(3)
    return add(x, mul(add(x, a), b)); //exp: s(7) = A + (A + A^4) * s(3)
}

int main() {
    while (scanf("%d%d%d", &n, &k, &m) == 3) {
        Mat res;
        memset(res.g, 0, sizeof(res.g));
        for (int i = 0; i < n; i++) 
            for (int j = 0; j < n; j++)
                scanf("%d", &res.g[i][j]);
        res = Cal(res, k);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n - 1; j++)
                printf("%d ", res.g[i][j]);
            printf("%d\n", res.g[i][n - 1]);
        }
    }
    return 0;
}