[HDU4471] Homework
问题描述
GS正苦恼于枯燥的数学作业。 这让他感到非常疲惫和不耐烦。为了能够出去溜达，好好享受生活，所以他让你来帮他完成数学作业。
在这项作业中, 你被要求解决以下问题：
给你一个递推公式：
和一些初始值
你需要求出\(f[n]\)。
多简单的问题啊。 等一等， 你看到了最后一段的几行。如下所述： 为了让这个问题难那么一点, 现在告诉你在一些特殊的位置(有q个位置, 被称为n1,n2,...,nq)， 通项公式会变成其他形式，即对于第k个特殊位置nk，通项公式变为
还是一个简单的问题，不是吗？
由于\(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项。不妨考虑构造矩阵来解决这个问题。来看一下这个递推式：
构造矩阵\(A_i\)为从\(i\)开始的前m项，那么
构造\(t\)行\(m\)列的转移矩阵\(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\n",++cnt,ans.a[1][1]);
}
return 0;
}
