【矩阵加速】[POJ3233]Matrix Power Series

题目描述

Description

Given a $n\times n$ matrix A and a positive integer k, find the sum $S = A + A^2 + A^3 + … + 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^9$) and m (m < $10^4$). Then follow n lines each containing n nonnegative integers below 32768, 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

题目分析

首先我们可以发现我们需要的就是每个矩阵的和那么我们用样例举例子

$$\begin{bmatrix}a_{11}&a_{12}&0&0 \\ a_{21}&a_{22}&0&0 \\ 0&0&0&0 \\ 0&0&0&0\end{bmatrix}\begin{bmatrix}a_{11}&a_{12}&1&0 \\ a_{21}&a_{22}&0&1 \\ 0&0&1&0 \\ 0&0&0&1\end{bmatrix}=\begin{bmatrix}a_{11}\times a_{12}+a_{12}\times a_{21}&...&a_{11}&a_{12} \\ ...&...&a_{21}&a_{22} \\ 0&0&0&0 \\ 0&0&0&0\end{bmatrix}$$
显然答案左上角为求得到的当前矩阵的多次方，右上角就是答案了，然后矩阵快速幂就好了

代码

#include <cstdio>
#include <iostream>
#include <cstring>
using namespace std;
const int MAXN = 53;
int MOD;
struct Matrix {
    int Ma[MAXN+10][MAXN+10];
    int n, m;
    void Clear(int u, int un, int um){
        n = un, m = um;
        for(int i=1;i<=un;i++)
            for(int j=1;j<=um;j++)
                Ma[i][j] = 0;
        if(u) for(int i=1;i<=un;i++)
            Ma[i][i] = 1;
    }
    Matrix operator* (const Matrix& ma) {
        Matrix ret ;
        ret.Clear(0, n, m);
        for(int i=1;i<=n;i++){
            for(int j=1;j<=ma.m;j++){
                for(int k=1;k<=ma.n;k++){
                    ret.Ma[i][j] += Ma[i][k] * ma.Ma[k][j];
                    ret.Ma[i][j] %= MOD;
                }
            }
        }
        return ret;
    }
    void Print(){
        for(int i=1;i<=n;i++){
            printf("%d", Ma[i][1]);
            for(int j=2;j<=m;j++)
                printf(" %d", Ma[i][j]);
            printf("\n");
        }
    }
};
Matrix Mpow(Matrix m, int p){
    Matrix ret;
    if(p == 0){
        ret.Clear(1, 2, 2);
        return ret;
    }else if(p == 1) return m;
    ret = Mpow(m, p/2);
    if(p%2 == 0) return ret * ret;
    return (ret * ret) * m;
}
int main(){
    int n, k;
    scanf("%d%d%d", &n, &k, &MOD);
    Matrix m, pm;
    m.Clear(0, n*2, n*2);
    pm.Clear(0, n*2, n*2);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            scanf("%d", &m.Ma[i][j]);
            pm.Ma[i][j] = m.Ma[i][j];
        }
        pm.Ma[i][i+n] = 1;
        pm.Ma[i+n][i+n] = 1;
    }
    //pm.Print();
    Matrix ans = Mpow(pm, k);
    Matrix tans = m * ans;
    for(int i=1;i<=n;i++){
        printf("%d", tans.Ma[i][1+n]);
        for(int j=2;j<=n;j++)
            printf(" %d", tans.Ma[i][j+n]);
        printf("\n");
    }

    return 0;
}