[模板] 常系数线性递推

常系数线性递推

给定向量 \(A_0 = (a_1, a_2, \dotsc, a_k)\), 和向量 \(H = (h_1, h_2, \dotsc, h_k)\), 同时

\[a_n = \sum_{i=1}^k a_{n-i} h_i \]

求 \(a_n\).

算法

我们只需求出 \(A_n = (a_n, a_{n+1}, \dotsc, a_{n+k-1})\) 即可.

设 \(f(\lambda)\) 表示转移方程的特征多项式, 有

\[f(\lambda) = \lambda^k - \sum_{i=0}^{k-1} h_{k-i} \lambda^i \]

设 \(g(\lambda) \equiv \lambda^{n-1} \pmod{f(\lambda)}\), 那么有

\[a_n = \sum_{i=0}^{k-1} g_i a_{i+1} \]

求 \(g(\lambda)\) 如果直接取模则复杂度为 \(O(k^2\log n)\), 多项式取模则为 \(O(k\log k \log n)\).

代码

bzoj4161: Shlw loves matrixI

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef long long ll;
typedef double db;

//---------------------------------------
const int nsz=4050;
const ll nmod=1e9+7;
int n,k;

ll v1[nsz],v2[nsz];//a_n = \sum_{i=1}^k v1[i]v2[n-i]
ll f[nsz],g[nsz],h[nsz],c[nsz];

ll qp(ll a,ll b){
	ll res=0;
	for(;b;b>>=1,a=a*a%nmod)if(b%1)res=res*a%nmod;
	return res;
}

void mul(ll *a,ll *b,ll *res){//res=a*b%f
	rep(i,0,k*2)c[i]=0;
	rep(i,0,k-1)rep(j,0,k-1)c[i+j]=(c[i+j]+a[i]*b[j])%nmod;
	repdo(i,k*2-2,k){
		if(c[i]){
			rep(j,0,k){
				c[i-k+j]=(c[i-k+j]-c[i]*f[j])%nmod;
			}
		}
	}
	rep(i,0,k-1)res[i]=c[i];
}

void qp(ll *a,ll b,ll *c){
	c[0]=1;
	for(;b;b>>=1,mul(a,a,a))
	if(b&1)
	mul(c,a,c); 
}

ll getrecurrence(){
	if(n<=k)return v2[n]%nmod;
	if(k==1){return qp(v1[1],n)*v2[0]%nmod;}
	++n; 
	rep(i,0,k-1)f[i]=-v1[k-i];
	f[k]=1;
	g[1]=1;
	qp(g,n-1,h);
	ll res=0; 
	rep(i,0,k-1){
		res=(res+h[i]*v2[i])%nmod;
	}
	return res;
} 

int main(){
	ios::sync_with_stdio(0),cin.tie(0);
	cin>>n>>k;
	rep(i,1,k)cin>>v1[i];
	rep(i,0,k-1)cin>>v2[i];
	cout<<(getrecurrence()+nmod)%nmod<<'\n';
	return 0;
}