矩阵快速幂模板题。（差别：这题目构建的矩阵中的的元素是数字， 下面题目是矩阵）
这个勇士明明超强却过分谨慎(矩阵快速幂， 模板)

正题
Matrix Power Series

题意：
给出$n,k,m$和一个$n\ast n$的矩阵$A$,求矩阵$S$,$S={\sum }_{i=1}^k(A^i)$.

$思路1$：
(1).用矩阵快速幂求$A^i$，再加起来。$TLE$了。
(2).根据公式找到递推式，再通过递推式构建新矩阵，用新矩阵做快速幂得到答案。
递推式：$S_i=A^i+S_{i-1}$
化成矩阵乘法会出现几种形式：(列举两种，其中$I$为单位矩阵)
1)$S_i=\left[\begin{array}{cc}{A}^i& S_{i-1}\end{array}\right]\ast \left[\begin{array}{c}I\\ I\end{array}\right]$

2)$S_i=\left[\begin{array}{c}{A}^{i-1}{S}_{i-1}\end{array}\right]\ast \left[\begin{array}{c}A\\ I\end{array}\right]$

在观察(1)式的时候，我们将右边也构造成$\left[\begin{array}{c}{A}^{i-1}{S}_{i-1}\end{array}\right]$的形式，$\left[\begin{array}{c}{A}^i{S}_i\end{array}\right]$，右边的$\left[\begin{array}{c}I\\ I\end{array}\right]$也要相应的改变，变成$\left[\begin{array}{cc}A& I\\ 0& I\end{array}\right]$。
合在一起就是$\left[\begin{array}{cc}{A}^{i+1}& S_i\end{array}\right]=\left[\begin{array}{cc}{A}^i& S_{i-1}\end{array}\right]\ast \left[\begin{array}{cc}A& I\\ 0& I\end{array}\right]$这样就和递推式的形式很像了。
同理(2)，也能化成$\left[\begin{array}{cc}{A}^i& S_i\end{array}\right]=\left[\begin{array}{cc}{A}^{i-1}& S_{i-1}\end{array}\right]\ast \left[\begin{array}{cc}A& A\\ 0& I\end{array}\right]$。

再用1)推出的公式可以变换成$\left[\begin{array}{cc}{A}^{i+1}& S_i\end{array}\right]=\left[\begin{array}{cc}A& 0\end{array}\right]\ast (\left[\begin{array}{cc}A& I\\ 0& I\end{array}\right]{)}^i$，这样一来我们就可以直接用一次矩阵快速幂和一个矩阵乘法就可以的出$S_i$矩阵

$Code$

// poj3233
#include<iostream>
#include<cstdio>
#include<vector>
#include<map>
#include<queue>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 35;
const ll INF = 0x7ffffffffffff;
const int inf = 0x7ffffff;
int n, k, m;
struct matrix {
    
    //矩阵，因为又要构建新矩阵，所以边要乘2
    ll s[2*N][2*N];
} I, ans, wor;
matrix mult(matrix V, matrix M) {
    
    //矩阵乘法
        matrix x; ll y, u = 2*n;
        memset(&x, 0, sizeof x);
        for(int i=1; i<=u; i++) {
    
    
            for(int j=1; j<=u; j++) {
    
    
                y = 0;
                for(int k=1; k<=2*n; k++) {
    
    
                    y = (y + V.s[i][k] * M.s[k][j]) % m;
                }
                x.s[i][j] = y % m;
            }
        }
        return x;
    }
matrix quick_mu(matrix A, ll u) {
    
       //矩阵快速幂（模板）
    if(u == 1) return A;
    matrix x = I;
    while(u) {
    
    
        if(u & 1) x = mult(x, A);
        A = mult(A, A);
        u >>= 1;
    }
    return x;
}
int read() {
    
    //快读（模板）
    int x = 0; int f = 1; char s = getchar();
    while(s < '0' || s > '9') {
    
    if(s == '-') f = -1; s = getchar();}
    while(s >= '0' && s <= '9') {
    
    x = (x << 3) + (x << 1) + s - 48; s = getchar();}
    return x * f;
}
int main() {
    
    
    for(int i=1; i<=2*N; i++) {
    
    
        I.s[i][i] = 1;//构建单位矩阵
    }
    n = read(), k = read(), m = read();
    for(int i=1; i<=n; i++) {
    
    
        for(int j=1; j<=n; j++) {
    
    
            ans.s[i][j] = read();//读入矩阵A
        }
    }
    wor = ans;
    for(int i=n+1; i<=2*n; i++) {
    
    //构建新矩阵
        wor.s[i][i] = 1;
        wor.s[i-n][i] = 1;
    }
    wor = quick_mu(wor, k);//用新矩阵进行快速幂
    ans = mult(ans, wor);//再做一次矩阵加法
    for(int i=1; i<=n; i++) {
    
    //输出答案
        for(int j=n+1; j<=2*n; j++) {
    
    
            printf("%lld ", ans.s[i][j]);
        }
        printf("\n");
    }
    return 0;
}

$思路2$:
(1).首先我们得会一般等比数列求和得公式。

//等比数列求和模板
ll qpow(ll a, ll b) {
    
    //快速幂：求a^b
    ll ans = 1; 
    while(b) {
    
    
        if(b & 1) ans = ans * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
}
ll faction(ll a, ll b) {
    
    //求：s_b;
    if(b == 0) return 1;
    if(b & 1) return faction(a, b/2) * (1+qpow(a, b/2+1)) % mod;
    else return (faction(a, (b-1)/2) * (1+qpow(a, (b-1)/2)) + qpow(a, b)) % mod;
}

(2).只要将上面两个函数中的元素都换成矩阵，乘法都换成矩阵的乘法，返回一个矩阵就行了。

matrix mqpow(matrix a, int x) {
    
    
    matrix res = I;//I为单位矩阵。
    while(x) {
    
    
        if(x & 1) res = mult(res, a);
        a = mult(a, a);
        x >>= 1;
    }
    return res;
}

matrix Y(matrix a, int x) {
    
    
    if(x == 0) return I;//返回单位矩阵。
    else if(x == 1) return a;
    else {
    
    
        if((x & 1) == 0) return mult(Y(a, x/2), fix(I, mqpow(a, x/2)));
        else return fix(mult(Y(a, (x-1)/2), fix(I, mqpow(a, (x-1)/2))), mqpow(a, x));
    }
}

$Code$

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
const int N = 35;
int n, k, mod;
typedef long long ll;
struct matrix {
    
    
    int m[N][N];
}A, ans, I;
void put(matrix s) {
    
    
    for(int i=1; i<=n; i++) {
    
    
        for(int j=1; j<=n; j++) {
    
    
            cout << s.m[i][j] << ' ';
        }
        cout << endl;
    }
}

matrix mult(matrix a, matrix b) {
    
    
    matrix x; ll y;
    memset(&x, 0, sizeof x);
    for(int i=1; i<=n; i++) {
    
    
        for(int j=1; j<=n; j++) {
    
    
            y = 0;
            for(int k=1; k<=n; k++) {
    
    
                y = (y + a.m[i][k]*b.m[k][j]) % mod;
            }
            
            x.m[i][j] = y % mod;
        }
    }
    return x;
}

matrix fix(matrix a, matrix b) {
    
    
    for(int i=1; i<=n; i++) {
    
    
        for(int j=1; j<=n; j++) {
    
    
            a.m[i][j] = (a.m[i][j]+b.m[i][j]) % mod;
        }
    }
    return a;
}
matrix mqpow(matrix a, int x) {
    
    
    matrix res = I;
    while(x) {
    
    
        if(x & 1) res = mult(res, a);
        a = mult(a, a);
        x >>= 1;
    }
    return res;
}

matrix Y(matrix a, int x) {
    
    
    if(x == 0) return I;
    else if(x == 1) return a;
    else {
    
    
        if((x & 1) == 0) return mult(Y(a, x/2), fix(I, mqpow(a, x/2)));
        else return fix(mult(Y(a, (x-1)/2), fix(I, mqpow(a, (x-1)/2))), mqpow(a, x));
    }
}

int main() {
    
    
    // freopen("in.txt", "r", stdin);
    cin >> n >> k >> mod;
    for(int i=1; i<=n; i++) I.m[i][i] = 1;
    for(int i=1; i<=n; i++) {
    
    
        for(int j=1; j<=n; j++) {
    
    
            cin >> A.m[i][j];
        }
    }
    put(Y(A, k));
    return 0;
}