BZOJ3231: [Sdoi2008]递归数列

BZOJ3231: [Sdoi2008]递归数列

Description

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

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

Input

由四行组成。

第一行是一个自然数k。

第二行包含k个自然数b₁, b₂,...,b_k。

第三行包含k个自然数c₁, c₂,...,c_k。

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

Output

仅包含一行：一个正整数，表示(a_m + a_m+1 + a_m+2 + ... + a_n) mod p的值。

Sample Input

2
1 1
1 1
2 10 1000003

Sample Output

142

HINT

对于100%的测试数据：
1<= k<=15
1 <= m <= n <= 1018

---

题解Here!

 

矩阵快速幂的沙茶题。

设$Ans(l,r)=\sum_{i=l}^ra_i$。

差一下分：$Ans(l,r)=Ans(1,r)-Ans(1,l-1)$

这种题只要构造出矩阵就万事大吉了。

我们很容易想到把$a_1,a_2,a_3,...,a_k$全部放到矩阵中。

但是求和怎么办？

没事，一并放到矩阵中。

设$sum(x)=\sum_{i=1}^xa_i$。

有这个式子：$$sum(x+1)=sum(x)+a_{x+1}=sum(x)+c_1\times a_{x}+c_2\times a_{x-1}+...+c_k\times a_{x-k+1}$$

所以我们构造出矩阵长这个样：$$\left[\begin{array}{}0&0&0&...&0&c_k&c_k\\1&0&0&...&0&c_{k-1}&c_{k-1}\\0&1&0&...&0&c_{k-2}&c_{k-2}\\0&0&1&...&0&c_{k-3}&c_{k-3}\\&&&......\\0&0&0&...&1&c_1&c_1\\0&0&0&...&0&0&1\end{array}\right]$$

最初的矩阵就是这样：$$\left[\begin{array}{}a_1&a_2&a_3&...&a_k&sum(k)\end{array}\right]$$

而$a_i=b_i,i\in [1,k]$。

然后就可以愉快地跑矩阵快速幂了。

附代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 20
using namespace std;
long long n,m,p,k;
long long b[MAXN],c[MAXN],sum[MAXN];
struct node{
	long long val[MAXN][MAXN];
	node(){
		for(int i=0;i<=19;i++)
		for(int j=0;j<=19;j++)
		val[i][j]=0;
	}
	friend node operator *(node x,node y){
		node ret;
		for(int i=1;i<=k+1;i++)
		for(int j=1;j<=k+1;j++){
			ret.val[i][j]=0;
			for(int l=1;l<=k+1;l++){
				ret.val[i][j]+=x.val[i][l]*y.val[l][j]%p;
				ret.val[i][j]%=p;
			}
		}
		return ret;
	}
	friend node operator ^(node x,long long w){
		node s;
		for(int i=1;i<=k+1;i++)s.val[i][i]=1;
		while(w){
			if(w&1)s=s*x;
			x=x*x;
			w>>=1;
		}
		return s;
	}
}a[3];
inline long long read(){
	long long date=0,w=1;char c=0;
	while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
	while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
	return date*w;
}
long long solve(long long x,int id){
	if(x<=k)return sum[x];
	node ans;
	for(int i=1;i<=k;i++)ans.val[1][i]=b[i];
	ans.val[1][k+1]=sum[k];
	ans=ans*(a[id]^(x-k));
	return ans.val[1][k+1]%p;
}
void work(){
	long long ans1=solve(n,1),ans2=solve(m-1,2);
	printf("%lld\n",(ans1-ans2+p)%p);
}
void init(){
	k=read();
	sum[0]=0;
	for(int i=1;i<=k;i++){
		b[i]=read();
		sum[i]=sum[i-1]+b[i];
	}
	for(int i=1;i<=k;i++)c[i]=read();
	m=read();n=read();p=read();
	a[1].val[k+1][k+1]=a[2].val[k+1][k+1]=1;
	for(int i=1;i<k;i++)a[1].val[i+1][i]=a[2].val[i+1][i]=1;
	for(int i=1;i<=k;i++)a[1].val[i][k]=a[1].val[i][k+1]=a[2].val[i][k]=a[2].val[i][k+1]=c[k-i+1];
}
int main(){
	init();
	work();
    return 0;
}