[bzoj 3231--SDOI2008]递归数列

一个由自然数组成的数列按下式定义：
对于i <= k：a[i]= b[i]
对于i > k: a[i] = c[1]*a[i-1] + c[2]*a[i-2] + …+ c[k]*a[i-k]
其中b[j]和 c[j] （1<=j<=k）是给定的自然数。写一个程序，给定自然数m <= n, 计算a[m] + a[m+1] + a[m+2] + … + a[n], 并输出它除以给定自然数p的余数的值。

这道题的算法显而易见是矩阵乘法与快速幂。只不过麻烦了一点，其余的地方还是比较基础的。讲一下矩阵构造吧，首先假设k=4，那么初始序列为a[1]~a[4],s[4]。s数组的意思挺好理解，举个例子，s[i]就表示a[1]+a[2]+…+a[i]。下图为相乘过程，图片中的第一个矩阵为半成品结果矩阵，由此可以得到构造半成品结果矩阵的规律。a[i][i+1]=1(1<=i,k>i),a[k][k-i+1]=a[k+1][k-i+1]=c[i](1<=i, k>=i),a[k+1][k+1]=1。为什么这样？图片已经告诉你了。得到构造半成品结果矩阵后，就很简单了。用两次快速幂，再乘初始序列，算出s[m-1]与s[n]，再相减，这题便解决了。当然有些情况需要特殊判断一下，所以这题才会麻烦。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdlib>
using namespace std;
struct node
{
    long long a[20][20];
    node()
    {
        memset(a,0,sizeof(a));
    }
};
long long n,p;
node chengfa1(node a,node b)
{
    node c;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            for(int k=1;k<=n;k++)
            {
                c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%p;
            }
        }
    }
    return c;
}
node chengfa2(node a,node b)
{
    node c;
    for(int i=1;i<=n;i++)
    {
        for(int k=1;k<=n;k++)
        {
            c.a[i][1]=(c.a[i][1]+a.a[i][k]*b.a[k][1])%p;
        }
    }
    return c;
}
node power(node a,long long x)
{
    node ans;
    for(int i=1;i<=n;i++)ans.a[i][i]=1;
    while(x>0)
    {
        if(x%2==1)ans=chengfa1(a,ans);
        a=chengfa1(a,a);
        x/=2;
    }
    return ans;
}
long long s[20];
int main()
{
    bool bk=true;
    long long l,r,x,s1,s2;
    node f,per,ans1,ans2;
    s[0]=0;
    scanf("%lld",&n);//n即为题目中的k
    for(int i=1;i<=n;i++)
    {
        scanf("%lld",&x);
        f.a[i][1]=x;f.a[n+1][1]+=x;
        s[i]=s[i-1]+x;
    }
    for(int i=1;i<=n;i++)//构造构造半成品结果矩阵
    {
        scanf("%lld",&x);
        if(i<n)per.a[i][i+1]=1;
        per.a[n][n-i+1]=per.a[n+1][n-i+1]=x;
    }
    per.a[n+1][n+1]=1;
    scanf("%lld%lld%lld",&l,&r,&p);
    l--;
    if(l<=n)s1=s[l];
    else
    {
        n++;
        bk=false;
        ans1=chengfa2(power(per,(l-n+1)),f);//其实应该l-n，但由于n++了，所以变成了l-(n-1)，下面关于r的一样的意思。
        s1=ans1.a[n][1];
    }
    if(bk==false)n--;
    if(r<=n)s2=s[r];
    else
    {
        n++;
        ans2=chengfa2(power(per,(r-n+1)),f);
        s2=ans2.a[n][1];
    }
    printf("%lld\n",((s2-s1)%p+p)%p);//一定要这样加p后再mod p，本题的一个大坑点
    return 0;
}