BZOJ3231 SDOI2008 递归数列 矩阵乘法

题意：设a_n=c₁a_n-1+c₂a_n-2+……+c_ka_n-k,s_n=a₁+a₂+……+a_n。给定a₁到a_k与c，求s_n-s_m-1

题解：以k=3为例\[\left( {\begin{array}{*{20}{c}}
{{a_{n - 2}}}&{{a_{n - 1}}}&{{a_n}}&{{s_n}}
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{a_1}}&{{a_2}}&{{a_3}}&{{s_3}}
\end{array}} \right){\left( {\begin{array}{*{20}{c}}
0&0&{{c_3}}&{{c_3}}\\
1&0&{{c_2}}&{{c_2}}\\
0&1&{{c_1}}&{{c_1}}\\
0&0&0&1
\end{array}} \right)^{n - 3}}\]

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;
#define ll long long

const int MAXN=16+2;
int K;
ll b[MAXN],c[MAXN],Sum[MAXN],Ans1,Ans2,P,M,N;
struct Matrix{
    int l,r;
    ll Mat[MAXN][MAXN];

    Matrix(){}
    Matrix(int _l,int _r):l(_l),r(_r){ memset(Mat,0,sizeof(Mat));}
    Matrix operator*(Matrix t){
        Matrix Ret=Matrix(l,t.r);
        for(int i=1;i<=l;i++)
        for(int j=1;j<=t.r;j++)
        for(int k=1;k<=r;k++)
            Ret.Mat[i][j]=(Ret.Mat[i][j]+Mat[i][k]*t.Mat[k][j]%P)%P;
        return Ret;
    }
    Matrix operator^(ll p){
        Matrix t=Matrix(l,r),Ret=Matrix(l,r);
        memcpy(t.Mat,Mat,sizeof(Mat));

        if(p&1) memcpy(Ret.Mat,Mat,sizeof(Mat));
        else
            for(int i=1;i<=l;i++) Ret.Mat[i][i]=1;

        while(p>>=1){
            t=t*t;
            if(p&1) Ret=Ret*t;
        }
        return Ret;
    }
}Int,Tran;

int main(){
    cin >> K;
    Int=Matrix(1,K+1);
    for(int i=1;i<=K;i++){
        scanf("%lld",b+i);
        Int.Mat[1][i]=b[i];
    }
    for(int i=1;i<=K;i++) scanf("%lld",c+i);
    scanf("%lld %lld %lld",&M,&N,&P);
    for(int i=1;i<=K;i++) Sum[i]=(Sum[i-1]+b[i])%P,Int.Mat[1][i]%=P,c[i]%=P;

    if(N<=K){
        cout << (Sum[N]-Sum[M-1]+P)%P << endl;
        return 0;
    }

    Tran=Matrix(K+1,K+1);
    for(int i=1;i<=K;i++) Tran.Mat[i][K]=Tran.Mat[i][K+1]=c[K-i+1];
    for(int i=2;i<=K;i++) Tran.Mat[i][i-1]=1;
    Tran.Mat[K+1][K+1]=1,Int.Mat[1][K+1]=Sum[K];

    if(M>K) Ans1=(Int*(Tran^(M-K-1))).Mat[1][K+1];
    else Ans1=Sum[M-1];
    Ans2=(Int*(Tran^(N-K))).Mat[1][K+1];
    cout << (Ans2-Ans1+P)%P << endl;

    return 0;
}