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

;
};
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;
}