[分治] 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\) 时间内完成矩阵乘法。

算法实现

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

转载于:https://www.cnblogs.com/lessmore/p/poj3233.html