【SDOI2008】【BZOJ3231】递归数列

Description

一个由自然数组成的数列按下式定义：
对于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的余数的值。
Input

由四行组成。
第一行是一个自然数k。
第二行包含k个自然数b1, b2,…,bk。
第三行包含k个自然数c1, c2,…,ck。
第四行包含三个自然数m, n, p。
Output

仅包含一行：一个正整数，表示(am + am+1 + am+2 + … + an) mod p的值。
Sample Input

2

1 1

1 1

2 10 1000003

Sample Output

142
HINT

对于100%的测试数据：

1<= k<=15

1 <= m <= n <= 1018

Source

矩乘+快速幂
把那一段和式看成两个前缀和相减就可以用矩乘了.

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define LL long long
#define MAXN 20
using namespace std;
int K;
LL b[20],c[20];
LL m,n,p;
LL ans;
struct matrix
{
    LL a[MAXN][MAXN];
    matrix()
    {
        memset(a,0,sizeof(a));
    }
    friend matrix operator *(matrix A,matrix B)
    {
        matrix ret;
        for (int i=1;i<=K+1;i++)
            for (int j=1;j<=K+1;j++)
                for (int k=1;k<=K+1;k++)
                    ret.a[i][j]=(ret.a[i][j]+A.a[i][k]*B.a[k][j])%p;
        return ret;
    }
    friend matrix operator ^(matrix x,LL k)
    {
        matrix ret;
        for (int i=1;i<=K+1;i++)    ret.a[i][i]=1;
        for (LL i=k;i;i>>=1,x=x*x)
        {
            if (i&1)    ret=ret*x;
        }
        return ret;
    }
}tmp1,tmp2;
LL solve(LL x)
{
    if (!x) return tmp2.a[1][K+1];
    matrix T=tmp1^x;
    T=tmp2*T;
    return T.a[1][K+1];
}
int main()
{
    scanf("%d",&K);
    for (int i=1;i<=K;i++)  scanf("%lld",&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++)  b[i]%=p,c[i]%=p,b[K+1]+=b[i],b[K+1]%=p;
    for (int i=1;i<=K;i++)  tmp1.a[i][1]=tmp1.a[i][K+1]=c[i];
    for (int i=2;i<=K;i++)  tmp1.a[i-1][i]=1;
    tmp1.a[K+1][K+1]=1;
    for (int i=1;i<=K;i++)  tmp2.a[1][i]=b[K-i+1];
    tmp2.a[1][K+1]=b[K+1];
    if (n<=K)
    {
        for (int i=m;i<=n;i++)  ans+=b[i],ans%=p;
        printf("%lld\n",ans);
        return 0;
    }
    ans=solve(n-K);
    if (m<=K)
        for (int i=1;i<m;i++)   ans-=b[i];
    else    ans-=solve(m-K-1);
    ans=(ans+p)%p;
    cout<<ans<<endl;
}