P2461 [SDOI2008]递归数列

题目描述

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

对于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的余数的值。

输入输出格式

输入格式：

输入文件spp.in由四行组成。

第一行是一个自然数k。

第二行包含k个自然数b1, b2,...,bk。

第三行包含k个自然数c1, c2,...,ck。

第四行包含三个自然数m, n, p。

输出格式：

输出文件spp.out仅包含一行：一个正整数，表示(am + am+1 + am+2 + ... + an) mod p的值。

输入输出样例

输入样例#1： 

2
1 1
1 1
2 10 1000003

输出样例#1： 

142

说明

对于100%的测试数据：

1<= k <=15

1 <= m <= n <= 10¹⁸

对于20%的测试数据：

1<= k <=15

1 <= m <= n <= 106

对于30%的测试数据：

k=1 1 <= m <= n <= 10¹⁸

对于所有测试数据：

0<= b1, b2,... bk, c1, c2,..., ck<=10⁹

1 <= p <= 10⁸

Solution：

　　本题矩阵快速幂。

　　求$sum_limits{i=m}^{ileq n}a_i$，可以转化为前缀和相减$s_n-s_{m-1}$。

　　那么我们需要快速求出$s_i$，我们发现$a_i$只与前$k$个$a$值有关，于是我们可以构建一个$(k+1)*(k+1)$的矩阵，存下前$k$个$a$值和当前的前缀和$s$。

　　转移矩阵的构造就补$1$并依次填好$c$值就好了。

代码：

/*Code by 520 -- 10.11*/
#include<bits/stdc++.h>
#define il inline
#define ll long long
#define RE register
#define For(i,a,b) for(RE int (i)=(a);(i)<=(b);(i)++)
#define Bor(i,a,b) for(RE int (i)=(b);(i)>=(a);(i)--)
using namespace std;
const int N=16;
struct matrix{
    int r,c;ll a[N][N];
    il void clr(){memset(a,0,sizeof(a));}
}ans,tp;
ll n,m,k,mod,b[N],c[N],s[N];

il matrix mul(matrix x,matrix y){
    matrix tp; tp.clr();
    tp.r=x.r,tp.c=y.c;
    For(i,0,x.r-1) For(j,0,y.c-1) For(k,0,x.c-1)
    tp.a[i][j]=(tp.a[i][j]+x.a[i][k]*y.a[k][j]%mod)%mod;
    return tp;
}

int main(){
    ios::sync_with_stdio(0);
    cin>>k;
    For(i,1,k) cin>>b[i],s[i]=(s[i-1]+b[i]);
    For(i,1,k) cin>>c[i];
    cin>>n>>m>>mod; ll tot=0;
    if(m<=k) cout<<(s[m]-s[n-1])%mod,exit(0);
    ans.r=1,ans.c=k+1; tp.r=tp.c=k+1; ans.clr(),tp.clr();
    For(i,0,k-1) ans.a[0][i]=b[i+1]%mod; ans.a[0][k]=s[k]%mod;
    For(i,1,k-1) tp.a[i][i-1]=1,tp.a[i][k-1]=tp.a[i][k]=c[k-i]%mod;
    tp.a[0][k-1]=tp.a[0][k]=c[k]%mod;tp.a[k][k]=1;
    if(n<=k) tot-=s[n-1]%mod;
    else {
        n-=k+1;
        while(n){
            if(n&1) ans=mul(ans,tp);
            n>>=1,tp=mul(tp,tp);
        }
        tot-=ans.a[0][k];
    }
    ans.r=1,ans.c=k+1; tp.r=tp.c=k+1; ans.clr(),tp.clr();
    For(i,0,k-1) ans.a[0][i]=b[i+1]%mod; ans.a[0][k]=s[k]%mod;
    For(i,1,k-1) tp.a[i][i-1]=1,tp.a[i][k-1]=tp.a[i][k]=c[k-i]%mod;
    tp.a[0][k-1]=tp.a[0][k]=c[k]%mod;tp.a[k][k]=1;
    m-=k;
    while(m){
        if(m&1) ans=mul(ans,tp);
        m>>=1,tp=mul(tp,tp);
    }
    tot=(tot+mod+ans.a[0][k])%mod;
    cout<<tot;
    return 0;    
}