常系数线性递推——多项式取模优化矩阵快速幂

题目：bzoj_4161

题目描述：

     给定数列 $h_n$前k项（下标从0开始），其后每一项满足$h_n=a_1\ast h_{n-1}+a_2\ast h_{n-2}+...+a_k\ast h_{n-k}$其中 a1,a2…ak 为给定数列。请计算 $h_n$，并将结果对 1000000007 取模输出。

题解：

     不妨将h的下标从1开始计，则所求值为$h_{n+1}$首先，最容易想出的办法是暴力，在暴力的基础上，我们可以引入矩阵快速幂，容易想到转移矩阵如下：
$$A=\left(\begin{array}{cccccc}{a}_1& a_2& a_3& \cdots & a_{k-1}& a_k\\ 1& 0& 0& \cdots & 0& 0\\ 0& 1& 0& \cdots & 0& 0\\ 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 1& 0\end{array}\right)$$
     而初始矩阵如下：
$$H_1=\left(\begin{array}{c}{h}_k\\ h_{k-1}\\ ⋮\\ h_2\\ h_1\end{array}\right)$$
     则转移过后的矩阵为
$$H_{n+1}=A^n\ast H_1=\left(\begin{array}{c}{h}_{k+n}\\ h_{k+n-1}\\ ⋮\\ h_{n+2}\\ h_{n+1}\end{array}\right)$$
     然而这样是会超时的。计算的瓶颈在于$A^n$的求法，因此需要用到矩阵的特征多项式取模来优化矩阵快速幂。将矩阵的乘法转成多项式的乘法。

前置知识：

1.特征值：

     设$A$是n阶矩阵，如果数$\lambda$和n维向量$x$使关系式$$Ax=\lambda x\text{\hspace{1em}(1)}$$成立，那么，这样的数$\lambda$称为矩阵$A$的 特征值

2.特征多项式：

     (1)式也可写成$$(A-\lambda E)x=0，$$这是n个未知数n个方程的齐次线性方程组，它有非零解的充分必要条件是系数行列式$$|A-\lambda E|=0,$$即
$$\left(\begin{array}{cccc}{a}_{11}-\lambda & a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}-\lambda & \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}-\lambda \end{array}\right)=0$$
上式是以$\lambda$为未知数的一元n次方程，成为矩阵$A$的***特征方程***，其左端$|A-\lambda E|$是$\lambda$的n次多项式，记作$f(\lambda )$，称为矩阵$A$的***特征多项式***。

3.Cayley-Hamilton定理：

     对$A$的特征多项式$f(\lambda )$，将矩阵$A$作为自变量带入，可得$$f(A)=0$$详细证明可以参考Wikipedia，这里是传送门。有一个很简单的理解方式就是直接代入，得到$f(A)=|A-AE|=0$，看起来还是挺显然的。

4.多项式取模：

     对于一个多项式$A(x)$，称其最高项的次数为这个多项式的度（degree），记作 $deg(A)$。设多项式$g(x)$为被除式，$deg(g)=n$，多项式$f(x)$为除式，$deg(f)=m（n>m)$。只要多项式$f(x)$不为常数0，则存在$$g(x)=f(x)\ast p(x)+r(x)$$，其中$deg(r)<m$。其中$p(x)$为商式，$r(x)$为余式。已知$f(x)$和$g(x)$来求$r(x)$的过程就叫做多项式取模。暴力取模可以直接模拟竖式除法。复杂度是$o(n^2)$（n与m同阶）。也可以用FFT来加速，达到log的复杂度，但是常数巨大，此题不适用。


     有了这四样东西后，我们就可以用来做题了。先看原来的转移矩阵的特征多项式，即
$$f(\lambda )=|A-\lambda E|=\left|\begin{array}{cccccc}{a}_1-\lambda & a_2& a_3& \cdots & a_{k-1}& a_k\\ 1& -\lambda & 0& \cdots & 0& 0\\ 0& 1& -\lambda & \cdots & 0& 0\\ 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 1& -\lambda \end{array}\right|$$

为了展开更方便，我们给等式乘上$(-1{)}^k$即
$$\begin{array}{rl}& (-1{)}^{k}f(\lambda )\\ =& (-1{)}^k|A-\lambda E|\\ =& (-1{)}^k\left|\begin{array}{cccccc}{a}_1-\lambda & a_2& a_3& \cdots & a_{k-1}& a_k\\ 1& -\lambda & 0& \cdots & 0& 0\\ 0& 1& -\lambda & \cdots & 0& 0\\ 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 1& -\lambda \end{array}\right|\\ =& \left|\begin{array}{cccccc}\lambda -a_1& -a_2& -a_3& \cdots & -a_{k-1}& -a_k\\ -1& \lambda & 0& \cdots & 0& 0\\ 0& -1& \lambda & \cdots & 0& 0\\ 0& 0& -1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & -1& \lambda \end{array}\right|\end{array}$$
将行列式按照第一行展开，可以得到
$$(-1{)}^{k}f(\lambda )={\lambda }^k-\sum _{i=1}^k{a}_i\ast x^{k-i}$$
不妨令该式为$g(\lambda )$，代入$A$，进行多项式取模，有如下等式
$$\begin{gathered} A^n=g(A)\ast p(A)+r(A) \\ deg(g)=k \end{gathered}$$则有$deg(r)<k$，不妨看作$deg(r)=k-1$，取模后得到r的系数$c_0,c_1,c_2\cdots c_{k-1}$。由于$g(A)=0$,得到$$A^n=r(A)=c_0+c_1\ast A^1+\cdots +c_{k-1}\ast A^{k-1}$$

则
$$H_{n+1}=A^n\ast H_1=\sum _{i=0}^{k-1}{c}_i\ast A_i\ast H_1=\sum _{i=0}^{k-1}{c}_i\ast H_{i+1}$$
即
$$H_{n+1}=c_0\ast \left(\begin{array}{c}{h}_k\\ ⋮\\ h_2\\ h_1\end{array}\right)+c_1\ast \left(\begin{array}{c}{h}_{k+1}\\ ⋮\\ h_3\\ h_2\end{array}\right)+\cdots +c_{k-1}\ast \left(\begin{array}{c}{h}_{2k}\\ ⋮\\ h_{k+1}\\ h_k\end{array}\right)$$
那么$h_{n+1}={\sum }_{i=0}^{k-1}{c}_i{h}_{i+1}$，即为答案。
求$c_n$数列的过程就是用$x^n$来模上$g(x)$所得的系数，这个过程可以用快速幂来实现。代码如下

#include<cstdio>
#include<algorithm>
using namespace std;
const int N=2005,mod=1000000007;
int n,k,a[N],h[N],c[N],b[N],tmp[N<<1];
void mul(int *a1,int *a2){
	for(int i=0;i<=2*k;i++)tmp[i]=0;
	for(int i=0;i<=k;i++){
		for(int j=0;j<=k;j++){
			tmp[i+j]=(tmp[i+j]+1ll*a1[i]*a2[j])%mod;
		}
	}
	for(int i=2*k;i>=k;i--){
		for(int j=0;j<k;j++){
			tmp[i-k+j]=(tmp[i-k+j]+1LL*tmp[i]*a[k-j])%mod;
		}
		tmp[i]=0;
	}
	for(int i=0;i<k;i++)a1[i]=tmp[i];
}
void pow(int *a1,int t,int *a2){
	for(;t;t>>=1){
		if(t&1)mul(a2,a1);
		mul(a1,a1);
	}
}
int main(){
#ifndef ONLINE_JUDGE
	freopen("1.in","r",stdin);
#endif
	scanf("%d%d",&n,&k);
	for(int i=1;i<=k;i++)scanf("%d",&a[i]);
	for(int i=1;i<=k;i++){
		scanf("%d",&h[i]);
		h[i]=(1ll*h[i]%mod+mod)%mod;
	}
	if(n<k)return printf("%d\n",h[n+1])*0;
	b[1]=1;c[0]=1;
	pow(b,n,c);
	int ans=0;
	for(int i=0;i<k;i++){
		ans=(ans+1LL*c[i]*h[i+1]%mod+mod)%mod;
	}
	printf("%d\n",ans);
	return 0;
}