一个由自然数组成的数列按下式定义：

 

对于i <= k：a_i = b_i

对于i > k: a_i = c₁a_i-1 + c₂a_i-2 + ... + c_ka_i-k

其中b_j和 c_j （1<=j<=k）是给定的自然数。写一个程序，给定自然数m <= n, 计算a_m + a_m+1 + a_m+2 + ... + a_n, 并输出它除以给定自然数p的余数的值。 1<= k<=15，1 <= m <= n <= 10^18

本题主要是构造矩阵，然后就直接矩阵连乘即可。

 

构造的上述矩阵的i大于k，当i小于k时就直接加起来就行。

 

这样a_m + a_m+1 + a_m+2 + ... + a_{n=S(n)-S(m-1)}

 

_代码：

#include <iostream>

#include <string.h>

#include <stdio.h>

using namespace std;

typedef long long LL;

const int MAXN=25;

struct Matrix

{

    LL m[MAXN][MAXN];

};

LL b[MAXN],c[MAXN];

LL K,M,N,P;

Matrix a,per;

void Init()

{

    int i,j;

    for(i=0;i<=K;i++)

        for(j=0;j<=K;j++)

            per.m[i][j]=(i==j);

}

Matrix multi(Matrix a,Matrix b)

{

    Matrix c;

    int i,j,k;

    for(i=0;i<=K;i++)

    {

        for(j=0;j<=K;j++)

        {

            c.m[i][j]=0;

            for(k=0;k<=K;k++)

                c.m[i][j]+=a.m[i][k]*b.m[k][j]%P;

            c.m[i][j]%=P;

        }

    }

    return c;

}

Matrix matrix_mod(LL n)

{

    Matrix ans=per,p=a;

    while(n)

    {

        if(n&1)

        {

            ans=multi(ans,p);

            n--;

        }

        n>>=1;

        p=multi(p,p);

    }

    return ans;

}

int main()

{

    int i,j;

    Init();

    LL ret1,ret2,S;

    Matrix ans;

    while(cin>>K)

    {

        ret1=ret2=0;

        for(i=0;i<K;i++)

           cin>>b[i];

        for(i=0;i<K;i++)

           cin>>c[i];

        cin>>M>>N>>P;

        for(i=0;i<K;i++)

        {

            b[i]%=P;

            c[i]%=P;

        }

        for(i=0;i<=K;i++)

        {

            for(j=0;j<=K;j++)

            {

                a.m[i][j]=0;

                if(i==0&&j==0) a.m[i][j]=1;

                if(i==0&&j>0)  a.m[i][j]=c[j-1];

                if(i==1&&j>0)  a.m[i][j]=c[j-1];

                if(i>1)        a.m[i][j]=(i==(j+1));

            }

        }

        S=0;

        for(i=0;i<K;i++)

        {

            S+=b[i];

            S%=P;

        }

        if(N<=K)

        {

            for(i=0;i<=N-1;i++)

            {

                ret1+=b[i];

                ret1%=P;

            }

        }

        else

        {

            ans=matrix_mod(N-K);

            ret1=ans.m[0][0]*S%P;

            for(i=1;i<=K;i++)

            {

                ret1+=ans.m[0][i]*b[K-i]%P;

                ret1%=P;

            }

        }

        if(M-1<=K)

        {

            if(M>=2)

            for(i=0;i<=M-2;i++)

            {

                ret2+=b[i];

                ret2%=P;

            }

            if(M<2) ret2=0;

        }

        else

        {

            ans=matrix_mod(M-K-1);

            ret2=ans.m[0][0]*S%P;

            for(i=1;i<=K;i++)

            {

                ret2+=ans.m[0][i]*b[K-i]%P;

                ret2%=P;

            }

        }

        cout<<((ret1-ret2)%P+P)%P<<endl;

    }

    return 0;

}