POJ 3233 Matrix Power Series (矩阵 + 二分)

Matrix67 博客上介绍的一道经典题目。

http://poj.org/problem?id=3233

题意 ： 求S = A + A^2 + A^3 + … + A^k 其中A是一个给定的n * n的矩阵， 对于m取模。n <= 30, m <= 10000, k <= 1e9

思路 ： 两次二分 ： 第一次二分计算calc(x) = A ^ x，第二次是S(1...x) = S(1...mid) +　S(mid + 1...x) = (calc(mid) + 1) * S(mid) + x是奇数则再加上calc(x)  注 ：mid = x / 2

CODE ：

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;
const int maxn = 35;
struct M {
    int a[maxn][maxn];
}bit[maxn], one, G;
int n, mod, K;

M mul(M A, M B) {
    M ans;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            ans.a[i][j] = 0;
            for (int k = 1; k <= n; k++) {
                ans.a[i][j] = (ans.a[i][j] + A.a[i][k] * B.a[k][j]) % mod;
            }
        }
    }
    return ans;
}
M add(M A, M B) {
    M ans;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            ans.a[i][j] = (A.a[i][j] + B.a[i][j]) % mod;
        }
    }
    return ans;
}
void init(M &a) {
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++)a.a[i][j] = 0;
}

M calc(int x) {
    M res = one;
    int c = 1;
    while (x) {
        if (x & 1) {
            res = mul(res, bit[c]);
        }
        c++; x >>= 1;
    }
    return res;
}

M solve(int x) { // calc(1...x)
    M res = one;
    if (!x) {
        return one;
    }else if (x == 1) {
        return G;
    }
    int mid = x >> 1;
    res = mul(add(one, calc(mid)), solve(mid));
    if (x & 1)res = add(res, calc(x));
    return res;
}
void Init() {
    bit[1] = G;
    for (int i = 2; i <= 33; i++) {
        bit[i] = mul(bit[i - 1], bit[i - 1]);
    }
}
void print(M a) {
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++)printf("%d ", a.a[i][j]);
        printf("\n");
    }
}

void BF() {
    M res; init(res);
    for (int i = 1; i <= K; i++) {
        res = add(res, calc(i));
    }
    printf("BF ::::\n");
    print(res);
}

int main() {
    scanf("%d%d%d", &n, &K, &mod);
    init(one); for (int i = 1; i <= n; i++)one.a[i][i] = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++){
            scanf("%d", &G.a[i][j]);
            G.a[i][j] %= mod;
        }
    }
    Init();
    print(solve(K));
    //BF();
    return 0;
}