[分治] POJ 3233 Matrix Power Series
题目大意
给出 \(n\times n\)矩阵 A 和正整数， 求出 $s= A + A^2 + A^3 + ... + A^k $。 取值范围， $n \leq 30, k \leq {10}^9, m < {10}^4 $.
解题思路
注意到 \(k\) 值很大，直接做矩阵乘法会超时，类似于二分快速幂算法，我们可以先求出 $s[k/2], A^{k/2} $，然后有下面递推式：
\[\begin{equation}
s[k]=
\begin{cases}
A^{k/2} * s[k/2] + s[k/2] &\mbox{if $k\bmod 2=0$ }\\
A^{k/2} * s[k/2] + s[k/2]+ A^k &\mbox{if $k \bmod 2 \neq 0 $ }
\end{cases}
\end{equation}
$$ \]
\begin{equation}
A^k = \begin{cases}
A^{k/2} * A^{k/2} &\mbox{if $ k\bmod 2=0 \(}\\
A^{k/2}*A^{k/2} * A &\mbox{if\) k \bmod 2 \neq 0 $ }
\end{cases}
\end{equation}
\[
这样就可以在 $\log k$ 时间内完成矩阵乘法。
## 算法实现
```c++
#include <stdio.h>
#include <string.h>
using namespace std;
int n, m;
int mat[30][30];
int rs[30][30]; // mat^1 + mat^2 + ... + mat^k
int r[30][30]; // mat^k
int tr[30][30]; // temp varible
void matmul(const int a[30][30], const int b[30][30], int d[30][30]) // mat multiple: d = a*b
{
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
d[i][j] = 0;
for (int k = 0; k < n; k++) {
d[i][j] += (a[i][k] * b[k][j]) % m;
}
d[i][j] %= m;
}
}
return;
}
void calc(int a) // after calc, r=mat^a, rs=mat^1 + ... + mat^a
{
if (a == 1) {
memcpy(rs, mat, sizeof(mat));
memcpy(r, mat, sizeof(mat));
return;
}
int b = a / 2;
calc(b);
// tr = rs*r + rs
matmul(rs, r, tr);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
tr[i][j] = (tr[i][j] + rs[i][j]) % m;
}
}
// rs <- tr
// tr = r*r
memcpy(rs, tr, sizeof(tr));
matmul(r, r, tr);
if (a % 2 == 0) {
memcpy(r, tr, sizeof(tr));
return;
}
else {
// r <- tr*mat
// rs += r
matmul(tr, mat, r);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
rs[i][j] = (rs[i][j] + r[i][j]) % m;
}
}
}
return;
}
int main()
{
int k;
scanf("%d %d %d", &n, &k, &m);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
scanf("%d", &mat[i][j]);
}
calc(k);
for (int i = 0; i < n; i++) {
printf("%d", rs[i][0]);
for (int j = 1; j < n; j++) {
printf(" %d", rs[i][j]);
}
printf("\n");
}
return 0;
}
```\]
