BZOJ3231(矩阵连乘，稍有点复杂)

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

对于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;
}