poj3233（等比矩阵求和）

poj3233

题意

给出一个 (n imes n) 的矩阵 (A) ，求 (A + A^2 + A^3 + ... + A^k) 。

分析

构造矩阵

[egin{bmatrix} A & E \ 0 & E \ end{bmatrix} ]

记为 (B) ，其中 (A) 为原矩阵，(E) 为 (n imes n) 的单位矩阵，(0) 为 (n imes n) 的零矩阵。
那么求 (B^{k+1}) ，
有

[egin{bmatrix} A^{k+1} & E+A+A^2+..+A^{k} \ 0 & E end{bmatrix} ]

可以发现右边减去一个单位矩阵就是我们所要求的答案。

code

#include<cstdio>
#include<cstring>
using namespace std;
const int N = 65;
int n, m, k;
struct Matrix {
    int mat[N][N];
    Matrix() { memset(mat, 0, sizeof mat); }
};
Matrix operator * (Matrix A, Matrix B) {
    Matrix C;
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            for(int k = 0; k < n; k++) {
                (C.mat[i][j] += A.mat[i][k] * B.mat[k][j]) %= m;
            }
        }
    }
    return C;
}
Matrix operator ^ (Matrix A, int x) {
    Matrix B;
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            if(i == j) B.mat[i][j] = 1;
        }
    }
    while(x) {
        if(x & 1) B = B * A;
        A = A * A;
        x >>= 1;
    }
    return B;
}
int main() {
    scanf("%d%d%d", &n, &k, &m);
    Matrix A;
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            scanf("%d", &A.mat[i][j]);
        }
    }
    // 右边的单位矩阵
    for(int i = 0; i < n; i++) {
        for(int j = n; j < 2 * n; j++) {
            if(j - i == n) A.mat[i][j] = 1;
        }
    }
    // 右下方的单位矩阵
    for(int i = n; i < 2 * n; i++) {
        for(int j = n; j < 2 * n; j++) {
            if(i == j) A.mat[i][j] = 1;
        }
    }
    // 下方的 0 矩阵，省略
    // ...

    n <<= 1;
    A = A ^ (k + 1);
    n >>= 1;

    // 将右边的矩阵减去一个单位矩阵
    for(int i = 0; i < n; i++) {
        for(int j = n; j < 2 * n; j ++) {
            if(j - i == n) A.mat[i][j] = (A.mat[i][j] - 1 + m) % m;
            printf("%d%c", A.mat[i][j], " 
"[j == 2 * n - 1]);
        }
    }
    return 0;
}