bzoj4161 Shlw loves matrix I 常系数齐次线性递推

Description

---

给定数列 {hn}前k项，其后每一项满足
hn = a1h(n-1) + a2h(n-2) + … + ak*h(n-k)
其中 a1,a2…ak 为给定数列。请计算 h(n)，并将结果对 1000000007 取模输出。

n <= 10⁹;k <= 2000; abs(hi)<=10⁹; abs(ai)<=10⁹

Solution

---

线性代数真是奥妙，推荐一篇讲得比较完整的blog

大概就是我们写出转移矩阵，求出特征方程，然后快速幂+多项式取模得到余柿，就可以用余柿算出答案了
如果套用FFT等操作多项式取模可以做到nlogn，直接暴力取模则是n^2的

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (register int i=st;i<=ed;++i)
#define drp(i,st,ed) for (register int i=st;i>=ed;--i)

typedef long long LL;
const int MOD=1000000007;
const int N=4005;

int n,m,k;
int mo[N],tmp[N],a[N],b[N],ans[N],h[N];

int read() {
	int x=0,v=1; char ch=getchar();
	for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):(v),ch=getchar());
	for (;ch<='9'&&ch>='0';x=x*10+ch-'0',ch=getchar());
	return x*v;
}

void upd(int &x,int v) {
	x+=v; (x>=MOD)?(x-=MOD):0;
}

void mul(int *a,int *b,int *c) {
	rep(i,0,k*2-2) tmp[i]=0;
	rep(i,0,k-1) rep(j,0,k-1) {
		upd(tmp[i+j],1LL*a[i]*b[j]%MOD);
	}
	drp(i,k*2-2,k) {
		drp(j,k-1,0) {
			upd(tmp[i+j-k],MOD-1LL*tmp[i]*mo[j]%MOD);
		}
	}
	rep(i,0,k-1) c[i]=tmp[i];
}

int main(void) {
	n=read(),k=read();
	rep(i,1,k) {
		a[i]=read();
		a[i]=(a[i]%MOD+MOD)%MOD;
	}
	rep(i,0,k-1) {
		h[i]=read();
		h[i]=(h[i]%MOD+MOD)%MOD;
	}
	if (n<k) return 0&printf("%d\n", h[n]);
	rep(i,0,k-1) mo[i]=(MOD-a[k-i])%MOD;
	if (k<=1) b[0]=-mo[0];
	else b[1]=1;
	mo[k]=1; ans[0]=1;
	for (n-=k-1;n;n>>=1) {
		if (n&1) mul(ans,b,ans);
		mul(b,b,b);
	}
	rep(i,k,k*2-1) {
		rep(j,1,k) upd(h[i],1LL*h[i-j]*a[j]%MOD);
	}
	int prt=0;
	rep(i,0,k-1) upd(prt,1LL*ans[i]*h[i+k-1]%MOD);
	prt=(prt%MOD+MOD)%MOD;
	printf("%d\n", prt);
	return 0;
}