【学习笔记】常系数线性齐次递推
应用(LGP4723)
求一个满足\(k\)阶线性齐次递推数列\(a_i\)的第\(n\)项；
及：\(a_n \ = \ \sum_{i=1}^{k} f_i \times a_{n-i}\) ;
$N \le 10^9 \ , \ K \le 32000 $ ；
矩阵
- 矩阵的秩\(rk(A)\)：
- \(rk(A) = 0 \ \Leftrightarrow \ A=0\) ;
- \(rk(A+B) \le rk(A) + rk(B)\) ;
- \(rk(AB) \le min(rk(A),rk(B))\) ;
- 行列式按行展开：
- \(det(A) \ = \sum_{j=1}^{n} a_{i,j}A_{i,j}\)
- \(A_{i,j}\)为行列式的余子式；
- 利用行列式的分拆性质可以证明；
特征多项式
- 对于\(n\)阶方阵\(A\)；
- 若存在列向量\(v\)和常数\(\lambda\) ，满足 \(\lambda \vec{v} = A \ \vec{v}\);
- 则称\(v\)为\(A\)的特征向量;
- 若\(k = rk(A)\)，则存在\(k\)个线性无关的\(v\)；
- 变换形式：\((\lambda I - A) \vec{v} = 0\) ;
- 有解的充要条件是\(det(\lambda I - A) = 0\);
- 所以可以得到一个关于 $ \lambda $ 的 $ n $ 阶方程 $ F(\lambda)=0 $ ， $ F(x) $ 称为 $ A $ 的特征多项式；
- Cayley-Hamilton定理
- 矩阵的特征多项式为矩阵的零化多项式：\(F(A) =0\) ；
- 证明：https://blog.csdn.net/qq_35649707/article/details/78688235 ；
线性递推
-
设递推数列 \(a_{0} - a_{k-1}\) 的初始列向量为\(G\)(下标从下到上)，\(A\)为\(k\)阶递推矩阵，实际上就是要求\([A^nG]_{k,1}\) ；
-
考虑如何求\(A^n\)，注意到\(F(A)=0\)；
-
设\(A^n = P(A) F(A)+ Q(A)\);
-
那么相当于\(A^n\)减去了很多个0矩阵，\(Q(A)\)和它意义相同，只需要算多项式\(Q(A)\)即可；
-
这部分可以快速幂+多项式取模完成，复杂度\(O(k log \ k log \ m)\)或者$k ^ 2 log \ n $；
-
考虑如何求\(A\)的特征多项式，利用定义$det(\lambda I - A) $的形式为:
\[\begin{pmatrix} \lambda - f_{1} & -f_2 & -f_3 & \cdots & -f_{k-1} & -f_{k} \\ 1 & \lambda & 0 & \cdots & 0 & 0 \\ 0 & 1 & \lambda & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & \lambda & 0 \\ 0 & 0 & 0 & \cdots & 1 & \lambda \\ \end{pmatrix}* \] -
按第一行展开，会剩下一个下三角行列式，值为对角线乘积；
-
得$F(\lambda) \ = \ \lambda^k - \sum_{i=1}^{k}f_i \lambda^{k-i} $ ;
-
\[g_{n} = [A^n G]_{k,1} = \sum_{i=0}^{k-1}[Q_{i} A^i \times G]_{k,1} =\sum_{i=0}^{k-1}Q_{i}[A^{i}G]_{k,1} = \sum_{i=0}^{k-1}Q_{i} a_{i} \]
-
总复杂度就是求\(Q\)的复杂度，另外如果前\(k\)项没有给出的话似乎可以利用生成函数求；
#include<bits/stdc++.h> #define mod 998244353 #define ll long long using namespace std; const int N=1000010; int n,k,m,f[N],h[N],a[N],b[N],ib[N]; char gc(){ static char*p1,*p2,s[1000000]; if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin); return(p1==p2)?EOF:*p1++; } int rd(){ int x=0,f=1;char c=gc(); while(c<'0'||c>'9'){if(c=='-')f=-1;c=gc();} while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+c-'0',c=gc();} return x*f; } int pw(int x,int y){ int re=1; if(y<0)y+=mod-1; while(y){ if(y&1)re=(ll)re*x%mod; y>>=1;x=(ll)x*x%mod; } return re; } void inc(int&x,int y){x+=y;if(x>=mod)x-=mod;} namespace poly{ const int G=3; int rev[N],L; void ntt(int*A,int len,int f){ for(L=0;(1<<L)<len;++L); for(int i=0;i<len;++i){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); if(i<rev[i])swap(A[i],A[rev[i]]); } for(int i=1;i<len;i<<=1){ int wn=pw(G,f*(mod-1)/(i<<1)); for(int j=0;j<len;j+=i<<1){ int w=1; for(int k=0;k<i;++k,w=(ll)w*wn%mod){ int x=A[j+k],y=(ll)w*A[j+k+i]%mod; A[j+k]=(x+y)%mod,A[j+k+i]=(x-y+mod)%mod; } } } if(!~f){ int iv=pw(len,mod-2); for(int i=0;i<len;++i)A[i]=(ll)A[i]*iv%mod; } } void cls(int*A,int l,int r){for(int i=l;i<r;++i)A[i]=0;} void cpy(int*A,int*B,int l){for(int i=0;i<l;++i)A[i]=B[i];} void inv(int*A,int*B,int l){ if(l==1){B[0]=pw(A[0],mod-2);return;} static int t[N]; int len=l<<1; inv(A,B,l>>1); cpy(t,A,l);cls(t,l,len); ntt(t,len,1);ntt(B,len,1); for(int i=0;i<len;++i)B[i]=(ll)B[i]*(2-(ll)t[i]*B[i]%mod+mod)%mod; ntt(B,len,-1);cls(B,l,len); } void pmod(int*A){ static int t[N]; int l=k+1,len=1;while(len<=(k<<1))len<<=1; cpy(t,A,(k<<1)+1); reverse(t,t+(k<<1)+1); cls(t,l,len); ntt(t,len,1); for(int i=0;i<len;++i)t[i]=(ll)t[i]*ib[i]%mod; ntt(t,len,-1); cls(t,l,len); reverse(t,t+l); ntt(t,len,1); for(int i=0;i<len;++i)t[i]=(ll)t[i]*b[i]%mod; ntt(t,len,-1); cls(t,l,len); for(int i=0;i<k;++i)A[i]=(A[i]-t[i]+mod)%mod; cls(A,k,len); } void pow(int*A,int n){ if(n==1){cls(A,0,k+1);A[1]=1;return;} pow(A,n>>1); int len=1;while(len<=(k<<1))len<<=1; ntt(A,len,1); for(int i=0;i<len;++i)A[i]=(ll)A[i]*A[i]%mod; ntt(A,len,-1); pmod(A); if(n&1){ for(int i=k;i;--i)A[i]=A[i-1];A[0]=0; pmod(A); } } } int main(){ // freopen("P4723.in","r",stdin); // freopen("P4723.out","w",stdout); n=rd();k=rd(); for(int i=1;i<=k;++i)f[i]=(mod+rd())%mod; for(int i=0;i<k;++i)h[i]=(mod+rd())%mod; for(int i=a[k]=b[k]=1;i<=k;++i)a[k-i]=b[k-i]=(mod-f[i])%mod; int len=1;while(len<=(k<<1))len<<=1; reverse(a,a+k+1); poly::inv(a,ib,len); poly::cls(ib,k+1,len); /*{ for(int i=0;i<len;++i)printf("%d ",b[i]);puts(""); for(int i=0;i<len;++i)printf("%d ",ib[i]);puts(""); }*/ poly::ntt(b,len,1); poly::ntt(ib,len,1); poly::pow(a,n); int ans=0; for(int i=0;i<k;++i)inc(ans,(ll)a[i]*h[i]%mod); printf("%d\n",ans); return 0; }
