问题描述
求一个满足 $K$ 阶齐次线性递推数列 $a_i$ 的第 $n$ 项，即：$a_n = sum_{i=1}^k f_i imes a_{n-i}$.
分析
首先写成矩阵快速幂
$$left( egin{bmatrix} f_1 &f_2 &f_3 &f_4 & cdots &f_{k-2} &f_{k-1} \ 1 &0 &0 &0 & cdots &0 &0 \ 0 &1 &0 &0 & cdots &0 &0\ cdots & cdots& cdots & cdots & cdots & cdots & cdots\ 0 &0 &0 &0 & cdots &1 &0 end{bmatrix} ight) ^n imes egin{bmatrix} a_{k-1} \ a_{k-2} \ cdots \ a_{1} \ a_{0}end{bmatrix} =egin{bmatrix} a_{n+k-1} \ a_{n+k-2} \ cdots \ a_{n+1} \ a_{n}end{bmatrix}$$
所以我们只需要算出 $M^N imes A$，然后取最后一个数即可。
使用矩阵快速幂，复杂度 $O(k^3 log_2n)$.
Carlay-Hamilton定理
设有 $k$ 个特征值的矩阵 $A$的特征多项式为 $f(lambda ) =prod_{i=1}^k(lambda_i - x)$，则有 $f(A) = 0$，$0$ 为零矩阵。
用这个定理来优化递推
由前面的矩阵快速幂，我们只要求出 $M^n$就可以了。
我们考虑 $M$ 的特征多项式 $f(x)$，这是一个 $k$ 次多项式。我们对 $M^n$ 做带余除法 $M^n = f(M) imes g(M) + R(M)$。
由于 $f(M) = 0$，所以 $M^n equiv R(M) (mod f(M))$，$R(M)$ 是一个次数不超过 $k-1$ 的多项式。
也就是说，我们只要求出 $M^n \% f(M)$就可以了
但是要怎么求呢？我们考虑快速幂的过程（就是倍增）
假设我们现在已知 $g(M)=M^{2^i} \% f(M)$，现在要求 $h(M)= M^{2^{i+1}} \% f(M)$。
一个直接的想法是令 $H(M)=g(M) imes g(M)$。但是这样做 $H(x)$ 的次数是 $2k-2$次的。
那么我们考虑原本的递推关系，$a_n=sumlimits_{i=1}^{k}a_{n-i}*f_i$
不难得到 $M^n=sumlimits_{i=1} ^{k} M^{n-i} imes f_{i}$
所以我们可以用这个式子将多余的系数都向前压一位。
这样我们可以得到一个 $O(k^2 log_2 n)$ 的做法。
那么有没有优化的余地呢？我们从倍增的过程入手，可以发现 $H(M) = g(M) imes g(M)$ 的过程可以用FFT/NTT加速至 $O(k log_2k)$。
现在只要解决压系数就可以了，把 $H(M)$ 模 $f(M)$ 即可。
我们的推导一直用到这个特征多项式 $f(x)$，如何求得呢？
根据定义， $f(x) = det(xI - M)$，得到
$$f(x) = |x I - M| = egin{bmatrix} x- a_1 & -a_2 & -a_3 & cdots & -a_{k - 2} & -a_{k - 1} & -a_k \ -1 & x & 0 & cdots & 0 & 0 & 0 \ 0 & -1 & x & cdots & 0 & 0 & 0 \ 0 & 0 & -1 & cdots & 0 & 0 & 0 \ vdots & vdots & vdots & ddots & vdots & vdots & vdots \ 0 & 0 & 0 & cdots & -1 & x & 0 \ 0 & 0& 0 & cdots & 0 & -1 & xend{bmatrix}$$
对第一行进行展开，得到
$$f(x) = (x - a_1)M_{11} + (-a_2)M_{12} + cdots + (-a_k)M_{1n} = x ^ k - a_1 x ^ {k - 1} - a_2x ^ {k - 2} - cdots - a_k$$
代码1：
$O(k log_2k log_2n)$的做法
思路其实就是去做一个类似快速幂的操作，然后把乘法改成多项式下的，取模也改成多项式下的
// luogu-judger-enable-o2 #include<cstdio> #include<algorithm> using namespace std; typedef long long ll; const ll mod=998244353; const int N=65536+10; int n;int k;int rv[20][N];ll rt[20][20];int Len;ll tr1[N];ll tr2[N];long long st[N];long long xs[N]; ll sg[N];ll a[N];ll res[N];ll irg[N];ll q[N];ll rf[N];int DL=-1;ll ans=0;ll ret[N]; inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;} inline void ntt(ll* a,int o,int len,int d)//ntt { for(int i=0;i<len;i++)if(i<rv[d][i])swap(a[i],a[rv[d][i]]); for(int k=1,j=1;k<len;k<<=1,j++) for(int s=0;s<len;s+=(k<<1)) for(int i=s,w=1;i<s+k;i++,w=w*rt[o][j]%mod) {ll a0=a[i];ll a1=a[i+k]*w%mod;a[i]=(a0+a1)%mod,a[i+k]=(a0+mod-a1)%mod;} if(o==1){ll inv=po(len,mod-2);for(int i=0;i<len;i++)(a[i]*=inv)%=mod;} } inline void poly_inv(ll* a,ll* b,int len)//求逆 { b[0]=po(a[0],mod-2); for(int k=1,j=0;k<=len;k<<=1,j++) { for(int i=0;i<k;i++)tr1[i]=a[i];for(int i=0;i<k;i++)tr2[i]=b[i]; ntt(tr1,0,k<<1,j);ntt(tr2,0,k<<1,j); for(int i=0;i<(k<<1);i++)b[i]=tr2[i]*(2+mod-tr1[i]*tr2[i]%mod)%mod; ntt(b,1,k<<1,j);for(int i=k;i<(k<<1);i++)b[i]=0; } } inline void poly_mod(ll* a)//取模 { int mi=(k<<1);while(a[--mi]==0);if(mi<k)return; for(int i=0;i<(Len<<1);i++)rf[i]=0;for(int i=0;i<=mi;i++)rf[i]=a[i]; reverse(rf,rf+mi+1);for(int i=mi-k+1;i<=mi;i++)rf[i]=0;ntt(rf,0,Len<<1,DL+1); for(int i=0;i<(Len<<1);i++)q[i]=(rf[i]*irg[i])%mod;ntt(q,1,(Len<<1),DL+1); for(int i=mi-k+1;i<=(Len<<1);i++)q[i]=0;reverse(q,q+mi-k+1);ntt(q,0,(Len<<1),DL+1); for(int i=0;i<(Len<<1);i++)(q[i]*=sg[i])%=mod;ntt(q,1,(Len<<1),DL+1); for(int i=0;i<k;i++)(a[i]+=mod-q[i])%=mod;for(int i=k;i<=mi;i++)a[i]=0; } int main() { for(int i=0;i<=15;i++) for(int j=0;j<(1<<(i+1));j++)rv[i][j]=(rv[i][j>>1]>>1)|((j&1)<<i); for(int t=2,j=1;j<=18;t<<=1,j++)rt[0][j]=po(3,(mod-1)/t); for(int t=2,j=1;j<=18;t<<=1,j++)rt[1][j]=po(332748118,(mod-1)/t); scanf("%d%d",&n,&k); for(Len=1;Len<=k;Len<<=1,DL++); //预处理 for(int i=1;i<=k;i++){scanf("%lld",&xs[i]);xs[i]=xs[i]<0?xs[i]+mod:xs[i];} for(int i=0;i<k;i++){scanf("%lld",&st[i]);st[i]=st[i]<0?st[i]+mod:st[i];} for(int i=1;i<=k;i++)sg[k-i]=mod-xs[i];sg[k]=1;for(int i=0;i<=k;i++)ret[i]=sg[i]; for(int i=0;i<=k;i++)rf[i]=sg[i];reverse(rf,rf+k+1);poly_inv(rf,irg,Len); for(int i=0;i<=k;i++)rf[i]=0;ntt(sg,0,Len<<1,DL+1);ntt(irg,0,Len<<1,DL+1);a[1]=1;res[0]=1; while(n)//快速幂 { if(n&1) { ntt(res,0,Len<<1,DL+1);ntt(a,0,Len<<1,DL+1); for(int i=0;i<(Len<<1);i++)(res[i]*=a[i])%=mod; ntt(res,1,Len<<1,DL+1);ntt(a,1,Len<<1,DL+1);poly_mod(res); }ntt(a,0,Len<<1,DL+1);for(int i=0;i<(Len<<1);i++)(a[i]*=a[i])%=mod; ntt(a,1,Len<<1,DL+1);poly_mod(a);n>>=1; } for(int i=0;i<k;i++)(ans+=res[i]*st[i])%=mod; printf("%lld",ans); return 0; }
代码2：
$O(k^2 log_2n)$的做法
#include<cstdio> #include<cstring> #include<cstdlib> #include<cctype> #include<cmath> #include<iostream> #include<algorithm> #include<vector> #include<set> #include<map> #include<queue> #include<stack> #include<cassert> typedef long long ll; typedef unsigned long long ull; using namespace std; const int P=1000000007; const int MAXN=4010; //2*k+10 int n,k,ans; int f[MAXN],h[MAXN]; struct Matrix{ //其实是多项式 int a[MAXN]; Matrix (){memset(a,0,sizeof a);} int& operator [] (const int &i) {return a[i];} int operator [] (const int &i) const {return a[i];} inline Matrix operator * (const Matrix &rhs) const { Matrix ret; for(int i=0;i<k;i++) for(int j=0;j<k;j++) (ret[i+j]+=1ll*a[i]*rhs[j]%P)%=P; for(int i=2*k-2;i>=k;ret[i--]=0) for(int j=1;j<=k;j++) //这里就是多项式取模优化的地方 (ret[i-j]+=1ll*ret[i]*f[j]%P)%=P; //可以认为是暴力向前压系数 return ret; } }res; Matrix ksm(Matrix a,int b) { Matrix ret; ret[0]=1; for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a; return ret; } int main() { scanf("%d%d",&n,&k); for(int i=1;i<=k;i++) scanf("%d",&f[i]),f[i]=f[i]>0?f[i]:f[i]+P; for(int i=0;i<k;i++) scanf("%d",&h[i]),h[i]=h[i]>0?h[i]:h[i]+P; if(n<k) printf("%d ",h[n]); res[1]=1;ans=0; res=ksm(res,n); for(int i=0;i<k;i++) ans=(ans+1ll*res[i]*h[i]%P)%P; printf("%d ",ans); }
参考链接：
1. https://www.luogu.org/problemnew/solution/P4723
2. https://www.luogu.org/blog/Zhang-RQ/chang-ji-shuo-ji-ci-xian-xing-di-tui-chu-tan