POJ 3233 Matrix Power Series

Matrix Power Series
Time Limit: 3000MS Memory Limit: 131072K
Total Submissions: 20659 Accepted: 8660
Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). 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

给一个矩阵A 求 S = A + A^2 + A^3 + … A^k
这个写得很详细
就是利用矩阵快速幂
然后对k进行二分

例如 k = 10的话
k = 10时:
S(10) = ( A^1+A^2+A^3+A^4+ A^5 ) + A^5 * ( A^1+A^2+A^3+A^4+A^5 ) = S(5) + A^5 * S(5)
然后计算k = 5：
S(5) = ( A^1+A^2 ) + A^3 + A^3 * ( A^1+A^2 ) = S(2) + A^3 + A^3 * S(2)
然后计算k = 2 :
k = 2 有 : S(2) = A^1 + A^2 = S(1) + A^1 * S(1)
思路就很清晰了 k要分奇偶讨论

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>
#include <vector>
#include <cmath>
#include <stack>
#include <string>
#include <map>
#include <set>
#define pi acos(-1)
#define LL long long
#define ULL unsigned long long
#define inf 0x3f3f3f3f
#define INF 1e18
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
using namespace std;
typedef pair<int, int> P;
const int maxn = 35;
const int maxm = 35;

int n, m, mod;
struct Matrix{
    int l, r;   // line | row
    LL mat[maxn][maxm];
    Matrix(){
        l = r = 0;
    }
    void initZero(){
        memset(mat, 0, sizeof(mat));
    }
    void initUnit(){
        for (int i = 0; i < maxn; i++)
            mat[i][i] = 1;
    }
    Matrix operator * (const Matrix &M) const{
        Matrix tmp;
        tmp.initZero();
        tmp.l = l; tmp.r = M.r;
        for (int i = 1; i <= l; i++){
            for (int j = 1; j <= r; j++){
                if (!mat[i][j]) continue;
                for (int k = 1; k <= M.r; k++){
//                    tmp.mat[i][k] += (mat[i][j] * M.mat[j][k]) % mod;// 多模一次会T....
                    tmp.mat[i][k] += mat[i][j] * M.mat[j][k]; 
                    tmp.mat[i][k] %= mod; 
                }
            }
        }
        return tmp;
    }
    Matrix operator + (const Matrix &M) const {
        Matrix tmp;
        tmp.initZero();
        tmp.l = l; tmp.r = r;
        for (int i = 1; i <= l; i++){
            for (int j = 1; j <= r; j++){
                tmp.mat[i][j] = (mat[i][j] + M.mat[i][j]) % mod;
            }
        }
        return tmp;
    }
    Matrix operator ^ (int k) const {
        Matrix tmp;
        Matrix M = *this;
        tmp.initZero();
        tmp.initUnit();
        tmp.l = tmp.r = l;
        while (k){
            if (k & 1)
                tmp = tmp * M;
            k >>= 1;
            M = M * M;
        }
        return tmp;
    }
    void print(){
        for (int i = 1; i <= l; i++){
            cout << mat[i][1];
            for (int j = 2; j <= r; j++)
                cout << " " << mat[i][j];
            cout << endl;
        }
    }
};
Matrix Sum(Matrix A, int k)
{
    if (k == 1) return A;
    else {
        Matrix a, b;
        a = Sum(A, k >> 1);
        if (k & 1){
            b = A ^ ((k>>1) + 1);
            return a + b + a * b;
        }
        else{
            b = A ^ (k>>1);
            return a + a * b;
        }
    }
}
int main(void)
{
    while (cin >> n >> m >> mod)
    {
        Matrix A;
        A.initZero();
        A.l = A.r = n;
        for (int i = 1; i <= n; i++){
            for (int j = 1; j <= n; j++){
                cin >> A.mat[i][j];
                A.mat[i][j] %= mod;
            }
        }
        Matrix ans = Sum(A, m);
        ans.print();
    }

    return 0;
}