bzoj#4161-Shlw loves matrixI【常系数线性齐次递推】

正题

题目链接:https://darkbzoj.tk/problem/4161

---

题目大意

给出序列\(a\)，和\(h\)的\(0\sim k-1\)项，满足

\[h_n=\sum_{i=1}^na_ih_{n-i} \]

求\(h_n\)。

\(1\leq n\leq 10^9,1\leq k\leq 2000\)

---

解题思路

显然\(k\leq 2000\)我们就不能矩阵乘法了，然后就去研究了一下常系数线性齐次递推，而且这题\(k^2\)能过所以不用\(O(n\log n)\)的多项式取模。

首先对于一个\(k\times k\)的矩阵\(A\)的特征多项式\(f(x)\)的定义是

\[f(x)=\sum_{i=0}^kdet(A-k\lambda)x^i \]

而我们知道对于这种递推式的特征多项式就是

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

然后有个Cayley-Hamilton定理就是说对于\(A\)的特征多项式\(f(x)\)满足\(f(A)=0\)。

如果按照这种特殊的特征多项式去理解的话挺好懂的，但是实际上的就不会证了。

那么我们如果要求一个转移矩阵\(A^n=A^n\mod f(A)\)，我们可以求出一个多项式\(r(x)=x^n\mod f(x)\)，然后直接把\(A\)带入这个多项式就好了。

这个东西要用到快速幂+多项式取模，\(O(k^2)\)的话实现起来很简单，考虑一个\(x^{k+i}(i\geq 0)\mod f(x)\)会发生什么，因为\(f(x)[k]=1\)所以\(\lfloor \frac{x^{k+i}}{f(x)}\rfloor=x^i\)相当于把\(x^i(f(x)-x^k)\)再丢回去就好了。

然后我们得到了一个多项式\(r(x)\)，考虑怎么计算答案，只需要计算\(A^{n-k}\vec{h}\)的最后一项，所以我们实际上求的\(r(x)=x^{n-k}\mod f(x)\)，记\(r(x)\)的\(i\)次系数为\(b_i\)，那么我们有

\[r(A)\vec{h}=\sum_{i=0}^{k-1}b_iA^{i}\vec{h} \]

所以我们算出\(h\)的\(k\sim 2k-1\)项然后就有

\[h_n=\sum_{i=0}^{k-1}r(x)[i]h_{k+i} \]

时间复杂度：\(O(k^2\log n)\)

值得一提的是如果转移矩阵\(A\)不是一个递推式的矩阵，我们也可以用拉格朗日插值插出它的特征多项式从而降低平均的复杂度。

---

code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=4100,P=1e9+7;
ll n,k,a[N],p[N],h[N],b[N],f[N],tmp[N];
void Mul(ll *a,ll *b,ll *c){
	memset(tmp,0,sizeof(tmp));
	for(ll i=0;i<k;i++)
		for(ll j=0;j<k;j++)
			(tmp[i+j]+=a[i]*b[j]%P)%=P;
	for(ll i=2*k-2;i>=k;i--)
		for(ll j=0;j<k;j++)
			(tmp[i+j-k]+=P-tmp[i]*p[j]%P)%=P;
	for(ll i=0;i<k;i++)c[i]=tmp[i];
	return;
}
signed main()
{
	scanf("%lld%lld",&n,&k);
	for(ll i=1;i<=k;i++)
		scanf("%lld",&a[i]),a[i]=(a[i]+P)%P,p[k-i]=P-a[i];
	for(ll i=0;i<k;i++)
		scanf("%lld",&h[i]),h[i]=(h[i]+P)%P;
	for(ll i=k;i<(k<<1);i++)
		for(ll j=1;j<=k;j++)
			(h[i]+=h[i-j]*a[j]%P)%=P;
	if(n<(k<<1))return printf("%lld\n",h[n])&0;
	b[1]=p[k]=f[0]=1;
	ll B=n-k,ans=0;
	while(B){
		if(B&1)Mul(f,b,f);
		Mul(b,b,b);B>>=1;
	}
	for(ll i=0;i<k;i++)
		(ans+=f[i]*h[i+k]%P)%=P;
	printf("%lld\n",ans);
	return 0;
}