POJ 3233 Matrix Power Series 【经典矩阵快速幂＋二分】

任意门：http://poj.org/problem?id=3233

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 28619		Accepted: 11646

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

Source

POJ Monthly--2007.06.03, Huang, Jinsong

题意概括：

给一个 N 维方阵 A ，求 A+A^2+A^3+ ... +A^k ，结果模 m；

解题思路：

矩阵快速幂解决矩阵幂运算（本质是二分优化）；

求前缀和二分：

比如，当k=6时，有：
    A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)
    应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A + A^2 + A^3，即可得到原问题的答案。

AC code：

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#define LL long long
using namespace std;
const int MAXN = 55;
int N, Mod;

struct mat
{
    int m[MAXN][MAXN];
}base;

mat mult(mat a, mat b)
{
    mat res;
    memset(res.m, 0, sizeof(res));
    for(int i = 0; i < N; i++)
    for(int j = 0; j < N; j++){
        if(a.m[i][j])
            for(int k = 0; k < N; k++)
                res.m[i][k] = (res.m[i][k] + a.m[i][j] * b.m[j][k])%Mod;
    }
    return res;
}

mat add(mat a, mat b)
{
    mat res;
    for(int i = 0; i < N; i++)
    for(int j = 0; j < N; j++)
        res.m[i][j] = (a.m[i][j] + b.m[i][j])%Mod;
    return res;
}

mat qpow(mat a, int k)
{
    mat ans;
    memset(ans.m, 0, sizeof(ans.m));
    for(int i = 0; i < N; i++) ans.m[i][i] = 1;

    while(k){
        if(k&1) ans = mult(ans, a);
        k>>=1;
        a=mult(a, a);
    }
    return ans;
}

mat solve(int K)
{
    if(K == 1) return base;
    mat res;
    memset(res.m, 0, sizeof(res.m));
    for(int i = 0; i < N; i++) res.m[i][i] = 1;

    res = add(res, qpow(base, K>>1));
    res = mult(res, solve(K>>1));
    if(K&1) res = add(res, qpow(base, K));

    return res;
}

int main()
{
    int K;
    mat ans;
    scanf("%d %d %d", &N, &K, &Mod);
    for(int i = 0; i < N; i++)
    for(int j = 0; j < N; j++){
        scanf("%d", &base.m[i][j]);
    }
    ans = solve(K);
    for(int i = 0; i < N; i++){
    for(int j = 0; j < N-1; j++)
            printf("%d ", ans.m[i][j]);
        printf("%d
", ans.m[i][N-1]);
    }
    return 0;
}