[HDU4471] Homework

问题描述

GS正苦恼于枯燥的数学作业。 这让他感到非常疲惫和不耐烦。为了能够出去溜达，好好享受生活，所以他让你来帮他完成数学作业。
在这项作业中, 你被要求解决以下问题：
给你一个递推公式：

[f[n]=sum_{i=1}^{t}c[i]f[n-i] ]

和一些初始值

[f[i],1 le i le m ]

你需要求出(f[n])。
多简单的问题啊。 等一等， 你看到了最后一段的几行。如下所述： 为了让这个问题难那么一点, 现在告诉你在一些特殊的位置(有q个位置, 被称为n1,n2,...,nq)， 通项公式会变成其他形式，即对于第k个特殊位置nk，通项公式变为

[f[n_k]=sum_{i=1}^{t_k}c_k f[n_k-i] ]

还是一个简单的问题，不是吗？
由于(f[n])的值可能会非常大，你只需要输出(f[n] mod 10^9+7)的值。

输入格式

有多组数据。

对每组数据, 第一行包括三个整数 n (m < n ≤ 10 9), m (1 ≤ m ≤ 100), q (0 ≤ q ≤ 100). 第二行包括 m 个整数(f[1],f[2],......,f[m])

接下来的一行包括几个整数,第一个为 t (t ≤ 100), 后面 t 个整数(c[1], c[2], . . . , c[t])。

接下来的q行描述了 q 个特殊位置的递推公式，每一行包括了几个整数 nk, tk (tk ≤ 100, tk < nk), ck[1], ck[2], . . . , ck[tk], 含义如题面所述。 特别地，如果i≠j，那么有ni ≠ nj 。

所有整数都是小于等于(10^9)的非负整数。

输入文件由EOF结尾。

保证所有数据合法。

输出格式

对每一组数据， 输出一行“Case X: Y”，其中 X 代表是第几组数据 (从1开始) ，Y 是该组数据的答案。

题目大意

给定一个数列的前m项以及递推公式，在某些特殊位置的递推公式单独给出，求数列的第n项。

来源

HDU 4417

解析

最简单的思路

直接从1到n计算每一项的值。但是n达到了(10^9)级别，直接循环并不能在给定的时间内解决问题。我们需要一个高效的递推算法。

矩阵乘法

不考虑特殊的项，对于每一次递推，我们需要知道前m项的值，以此得到该项的值。也就是说，每次需要保留前m项。不妨考虑构造矩阵来解决这个问题。来看一下这个递推式：

[f[n]=sum_{i=1}^{t}c[i]f[n-i] ]

构造矩阵(A_i)为从(i)开始的前m项，那么

[A_i=left [ egin{matrix} f[i]&f[i-1]&...&f[i-m+1] end{matrix} ight ] ]

构造(t)行(m)列的转移矩阵(S)，利用矩阵乘法的规则有如下构造方案：

[S=left [ egin{matrix} c[1]&1&0&0&0&...&0\c[2]&0&1&0&0&...&0\.\.\.\c[t-1]&0&0&0&0&...&1\c[t]&0&0&0&0&...&0 end{matrix} ight ] ]

手动模拟一下就可以明白这样构造的原理。所以我们有

[A_{i+1}=A_i*S ]

然后呢？

矩阵快速幂

用矩阵快速幂就可以在(log)级别的时间复杂度内解决了。但是还有特殊位置的问题没有解决。

解决问题

不妨考虑将每一个特殊位置拎出来单独处理。对于两个特殊位置(p_1)和(p_2)之间用矩阵快速幂处理，计算(p_2)时带入公式处理。

但是直接矩阵快速幂还是会超时。优化一下，把每(2^i)个转移矩阵的乘积预处理出来。

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#define int long long
#define N 105
using namespace std;
const int mod=1000000007;
struct Matrix{
	int n,m,a[N][N];
}ans,s[N];
struct event{
	int n,t,c[N];
}a[N];
int n,m,p,t,size,i,j,f[N],c[N],cnt;
int read()
{
	char c=getchar();
	int w=0;
	while(c<'0'||c>'9') c=getchar();
	while(c<='9'&&c>='0'){
		w=w*10+c-'0';
		c=getchar();
	}
	return w;
}
int my_comp(const event &x,const event &y)
{
	return x.n<y.n;
}
Matrix mult(Matrix a,Matrix b)
{
	Matrix c;
	c.n=a.n;c.m=b.m;
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++) c.a[i][j]=0;
	}
	for(int i=1;i<=a.n;i++){
		for(int j=1;j<=b.m;j++){
			for(int k=1;k<=a.m;k++) c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
		}
	}
	return c;
}
Matrix poww(Matrix a,int b)
{
	for(int i=0;(1<<i)<=b;i++){
		if(b&(1<<i)) a=mult(a,s[i]);
	}
	return a;
}
signed main()
{
	while(scanf("%lld%lld%lld",&n,&m,&p)!=EOF){
		memset(s[0].a,0,sizeof(s[0].a));
		memset(ans.a,0,sizeof(ans.a));
		for(i=1;i<=m;i++) f[i]=read();
		size=t=read();
		for(i=1;i<=t;i++) c[i]=read();
		for(i=1;i<=p;i++){
			a[i].n=read();a[i].t=read();
			if(a[i].n<=n) size=max(size,a[i].t);
			for(j=1;j<=a[i].t;j++) a[i].c[j]=read();
		}
		sort(a+1,a+p+1,my_comp);
		s[0].n=s[0].m=size;
		for(i=1;i<size;i++) s[0].a[i][i+1]=1;
		for(i=1;i<=t;i++) s[0].a[i][1]=c[i];
		for(i=1;(1<<i)<=n;i++) s[i]=mult(s[i-1],s[i-1]);
		ans.n=1;ans.m=size;
		for(i=1,j=m;i<=size;i++){
			ans.a[1][i]=f[j];
			j--;
			if(!j) break;
		}
		int pos=m;
		for(i=1;i<=p;i++){
			if(a[i].n<=pos||a[i].n>n) continue;
			ans=poww(ans,a[i].n-pos-1);
			pos=a[i].n;
			int tmp=0;
			for(j=1;j<=a[i].t;j++) tmp=(tmp+a[i].c[j]*ans.a[1][j])%mod;
			for(j=ans.m;j>1;j--) ans.a[1][j]=ans.a[1][j-1];
			ans.a[1][1]=tmp;
		}
		ans=poww(ans,n-pos);
		printf("Case %lld: %lld
",++cnt,ans.a[1][1]);
	}
	return 0;
}