[分治] 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; } ```\]