C++题解：Matrix Power Series ——矩阵套矩阵的矩阵加速

Matrix Power Series

r时间限制: 1 Sec 内存限制: 512 MB
题目描述
给定矩阵A，求矩阵S=A^1+A2+……+A^k，输出矩阵，S矩阵中每个元都要模m。

数据范围： n (n ≤ 30) , k (k ≤ 109) ，m (m < 104)

输入
输入三个正整数n，k，m

输出
输出矩阵S mod m

样例输入

2 2 4
0 1
1 1
样例输出
1 2
2 3

这道题不多说，可以得出加速矩阵(E为单位矩阵，也就是形为$\begin{bmatrix}1&0&...&0\\0&1&...&0\\... &...&...&...\\0&0& ...&1\end{bmatrix}$的矩阵，任何矩阵乘以这个单位矩阵还是原矩阵）：
$\begin{bmatrix} A &E \\ 0 & E \end{bmatrix}$*$\begin{bmatrix} A &E \\ 0 & E \end{bmatrix}$=$\begin{bmatrix} A^{2} &E+A \\ 0 & E \end{bmatrix}$
所以根据题目的要求，答案便是$\begin{bmatrix} A &E \\ 0 & E \end{bmatrix}^{k+1}$的(1,2)
主要难点是矩阵套矩阵，详见代码：

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
      
#define N 35
#define P 5
#define LL long long
     
LL fuck,k,mod;
 
struct M {
    int n,m,c[N][N];
    M() { 
        n=m=fuck; 
        memset(c,0,sizeof(c)); 
    } 
    M operator * (const M& a) {
        M r;
        r.n=n;r.m=a.m;
        for(int i=1;i<=r.n;i++)
            for(int j=1;j<=r.m;j++)
                for(int k=1;k<=m;k++)
                    r.c[i][j]= ( r.c[i][j] + (c[i][k] * a.c[k][j] ) % mod) % mod;
        return r;
    }
    M operator + (const M& a) {
        M r;
        for(int i=1;i<=r.n;i++)
            for(int j=1;j<=r.m;j++)
                r.c[i][j]=(c[i][j]+a.c[i][j]) %mod;
        return r;   
    }
    M operator - (const M& a) {
        M r;
        for(int i=1;i<=r.n;i++)
            for(int j=1;j<=r.m;j++)
                r.c[i][j]=r.c[i][j]+(mod+c[i][j]-a.c[i][j]) %mod;
        return r;
    }
    void _read() {
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
                scanf("%lld",&c[i][j]);
    }
    void pre() {
        n=m=fuck;
        for(int i=1;i<=fuck;i++)
            c[i][i]=1;
    }
    void _print() {
        for(int i=1;i<=n;i++) {
            for(int j=1;j<=m;j++) {
                if(j!=1) cout<<" ";
                cout<<c[i][j];
            }
            if(i!=n) puts("");
        }
        puts("");
    }
}fuckk;
 
struct Matrix {
    LL n,m;
    M c[P][P];
    Matrix() { 
        m=2,n=2;
        memset(c,0,sizeof(c)); 
    };
    Matrix operator * (const Matrix& a) {
        Matrix r;
        r.n=n;r.m=a.m;
        for(int i=1;i<=r.n;i++)
            for(int j=1;j<=r.m;j++)
                for(int k=1;k<=m;k++)
                    r.c[i][j]=r.c[i][j] + (c[i][k] * a.c[k][j] );
        return r;
    }
     
    Matrix pow(Matrix a, LL indexx) {
        Matrix sum;sum.n=sum.m=2;
        sum.c[1][1].pre();
        sum.c[2][2].pre();
        //a.c[1][2]._print();
        while(indexx>0) {
			if(indexx&1) sum=sum*a;
			/*sum.c[1][1]._print();
            sum.c[1][2]._print();
            sum.c[2][1]._print();
            sum.c[2][2]._print();*/
            a=a*a;
            //a.c[1][1]._print();
            indexx/=2;
        }
        return sum;
    }
    void sub() {
        c[1][2]=c[1][2]-fuckk;
    }
     
}ans;
  
int main() {
    cin>>fuck>>k>>mod;
    M a,b;
    a._read(); 
    b.pre();
    fuckk=b;
    ans.c[1][1]=a;
    ans.c[1][2]=ans.c[2][2]=b;
    //ans.test(ans);
    ans=ans.pow(ans,k+1);
    //ans.c[1][2]._print();
    ans.sub();
    ans.c[1][2]._print();
}