常系数线性齐次递推
给定数列前k项\(h_0…h_{k-1}\)，其后的项满足：\(h_i=\sum_{i=1}^kh_{i-j}a_i\)，其中\(a_1…a_k\)是给定的系数，求\(h_n\)
数据范围小的时候：
做法一：暴力\(O(nk)\)的DP
做法二：矩阵快速幂.
记\(H_i=\begin{bmatrix}h_i&h_{i+1}&…&h_{i+k-1}\end{bmatrix}\). 则\(h_n\)是\(H_{n-k+1}\)的最后一项。
\(H_{n-k+1}=H_0M^{n-k+1}\)
其中\(M\)是转移矩阵，如当\(k=4\)时是这么填的：
\[
M=\begin{bmatrix}
0&0&0&a_4\\
1&0&0&a_3\\
0&1&0&a_2\\
0&0&1&a_1
\end{bmatrix}
\]
时间复杂度\(O(k^3lg n)\)
数据范围大一些的时候：
\(k\leq2000,n\leq10^9\). 这时候矩阵快速幂也做不了了
还是拿\(k=4\)时举例，\(M\)的特征多项式\(f(\lambda)\)为：
\[
f(\lambda)=det(\lambda I-M)=\begin{bmatrix}
\lambda&0&0&0\\
0&\lambda&0&0\\
0&0&\lambda&0\\
0&0&0&\lambda
\end{bmatrix}
-\begin{bmatrix}
0&0&0&a_4\\
1&0&0&a_3\\
0&1&0&a_2\\
0&0&1&a_1
\end{bmatrix}=\begin{bmatrix}
\lambda&0&0&-a_4\\
-1&\lambda&0&-a_3\\
0&-1&\lambda&-a_2\\
0&0&-1&\lambda-a_1
\end{bmatrix}
\]
用行列式的性质，将\(f(\lambda)\)按最后一列拉普拉斯展开，得到如下，其中\((-1)^{i+j}f(x)_{i,j}\)即行列式定义里的代数余子式：
\[
\begin{aligned}
f(\lambda)&=\sum_{i=1}^ka_{k-i+1}(-1)^{i+j}f(\lambda)_{i,j} &取j=k（按最后一列展开）\\
&=\sum_{i=1}^ka_{k-i+1}(-1)^{i+k}f(\lambda)_{i,k}
\end{aligned}
\]
化简得到如下式子（也可以按\(k=4\)带进去看看规律）
\[
f(\lambda)=\lambda^k-\sum_{i=1}^ka_i\lambda^{k-i}
\]
现在明确一个定义，\(f(x)\)这个函数的自变量\(x\)可以是实数，也可以是矩阵等等。这个函数仅仅是表示如何将自变量组合起来。表达的意思也会多样化，比如多项式、矩阵的多项式…下文会随时切换自变量的种类，但是函数的本质不变。
\(\lambda\)是\(M\)的特征值，是一个数。但是根据Cayley-Hamilton定理，如果把\(\lambda\)替换成\(M\)代入得到\(f(M)=M^k-\sum_{i=1}^ka_iM^{k-i}\)，结果为一个零矩阵，即\(M^k-\sum_{i=1}^ka_iM^{k-i}=0\)
我们想要求\(M\)的\(n\)次方（这里的\(n\)只是代表\(M\)的\(n\)次方，题目中\(n\)应该用\(n-k+1\)替代），然而\(M^n\)直接快速幂求不现实，复杂度为\(O(k^3lg n)\).
首先退一步考虑，要求一个数字的n次方\(x^n\)，如果我们把\(x^n\)对\(f(x)\)取模会发生什么？
根据多项式取模的定义，\(x^n modf(x)=f(x)g(x)+r(x)\)，其中\(g(x)\)和\(r(x)\)是两个多项式（自变量为矩阵的多项式）.
将\(x\)看成\(M\)，那么\(f(M)\)为0.
故\(M^nmodf(M)=r(M)\)，且\(M^n=M^nmodf(M)\)，那么\(M_n=r(M)\)这个多项式
根据多项式取模的特性，\(r(x)\)的阶数严格小于模数\(f(x)\)的阶数\(k\)(最高次项是\(x^k\)). 那么\(r(x)\)所包含的\(M\)的指数一定小于\(k\)，到达了可以计算的范围。
要求\(M^n\)，就只需要求\(M^n mod f(M)\)的多项式\(r(M)\)。如果两个多项式\(A(x)\)和\(B(x)\)对模数取模分别得到\(C(x)\)和\(D(x)\)，那么多项式\(A(x)B(x)\)对模数取模结果就是\(C(x)D(x)\)。
那么就可以用快速幂来求解\(M^nmodf(M)\)的结果了，也就是求出了\(r(x)\)的各项系数（记为\(c_i\)）。实际计算中，表面上是在计算\(M^n\)，实际上计算的是\(M^nmodf(M)\)的结果。
至此求出\(r(x)=\sum\limits_{i=0}^{k-1}c_ix^i\). 将它看成矩阵的多项式代入\(M\)，得\(r(M)=\sum\limits_{i=0}^{k-1}c_iM^i\)
所以\(M^n=\sum\limits_{i=0}^{k-1}c_iM^i\)
把\(n\)替换成题目所需要的\(n-k+1\)，最终答案\(h_n\)为\(H_0M^{n-k+1}\)的最后一项。
\[
H_0M^{n-k+1}=H_0\sum_{i=0}^{k-1}c_iM^i=\sum_{i=0}^{k-1}c_iH_0M_i=\sum_{i=0}^{k-1}c_iH_i
\]
额那么要求的是\(H_0M^{n-k+1}\)的最后一项。记\(last(H_i)=h_{k+i}\) ，那么
\[
h_n=last(H_0M^{n-k+1})=\sum_{i=0}^{k-1}c_ilast(H_i)=\sum_{i=0}^{k-1}c_ih_{i+k}
\]
发现\(i+k\in[k,2k-1]\)，所以暴力算出\(h_k…h_{2k-1}\)，代入求解得到\(h_n\)，至此全部求完。
分析复杂度：多项式乘法此处用暴力\(n^2\)算会比FFT快，耗时最多的集快速幂求\(r(x)\) ，复杂度为\(O(n^2lgn)\)。
#include <cstdio>
using namespace std;
const int K=4005,mod=1e9+7;
int n,k;
int a[K],h[K];
int b[K],c[K],t[K],mo[K];
inline void add(int &x,int y){
x+=y;
if(x>=mod) x-=mod;
}
void mul(int *x,int *y,int *z){
for(int i=0;i<=2*k-2;i++) t[i]=0;
for(int i=0;i<k;i++)
for(int j=0;j<k;j++)
add(t[i+j],1LL*x[i]*y[j]%mod);
for(int i=2*k-2;i>=k;i--){
for(int j=k-1;j>=0;j--)
add(t[i-k+j],mod-1LL*t[i]*mo[j]%mod);
t[i]=0;
}
for(int i=0;i<k;i++) z[i]=t[i];
}
void ksm(int y){
for(;y;mul(b,b,b),y>>=1)
if(y&1)
mul(c,b,c);
}
int main(){
freopen("input.in","r",stdin);
scanf("%d%d",&n,&k); n++;
for(int i=1;i<=k;i++){
scanf("%d",&a[i]);
if(a[i]<0) a[i]+=mod;
}
for(int i=1;i<=k;i++){
scanf("%d",&h[i]);
if(h[i]<0) h[i]+=mod;
}
mo[k]=1;
for(int i=1;i<=k;i++) mo[k-i]=mod-a[i];
if(n<=k){printf("%d\n",h[n]);return 0;}
b[1]=1; c[0]=1;
ksm(n-k);
for(int i=k+1;i<=2*k;i++)
for(int j=1;j<=k;j++)
h[i]+=1LL*a[j]*h[i-j]%mod,h[i]-=h[i]>=mod?mod:0;
int ans=0;
for(int i=0;i<k;i++)
ans+=1LL*c[i]*h[i+k]%mod,ans-=ans>=mod?mod:0;
printf("%d\n",ans);
return 0;
}