BZOJ3231(矩陣連乘，稍有點復雜)

題意：

一個由自然數組成的數列按下式定義：
 
對於i <= k：ai = bi
對於i > k: ai = c1ai-1 + c2ai-2 + ... + ckai-k
其中bj和 cj （1<=j<=k）是給定的自然數。寫一個程序，給定自然數m <= n, 計算am + am+1 + am+2 + ... + an, 並輸出它除以給定自然數p的余數的值。1<= k<=15，1 <= m <= n <= 10^18

本題主要是構造矩陣，然後就直接矩陣連乘即可。

構造的上述矩陣的i大於k，當i小於k時就直接加起來就行。
 
這樣am + am+1 + am+2 + ... + an=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;
}