常系数齐次线性递推快速算法学习笔记
今天集训被线代狠狠的虐了一发。
不过还有一点收获的，比如这个。
数列 \(f\) 满足 \(f_n=\sum\limits_{i=1}^ka_if_{n-i}(n\ge k)\)，其中 \(a_1\dots a_k,f_0\dots f_{k-1}\) 均给出。求 \(f_n\)。
\(n\le 10^9,k\le 30000\)。
先要弄懂一些基本定义，
矩阵，行列式，高斯消元这些基本的东西就自己看别的东西去吧，我也不知道有什么好的资料，不妨去洛谷搜模板看题解，那里面的都不错。
然后讲一下特征值，特征多项式和 Hamilton-Cayley 定理，就能做这题了。
对于 \(n\times n\) 的矩阵 \(A\)，如果存在数 \(\lambda\) 和非零列向量 \(x\) 满足 \(Ax=\lambda x\)，那么 \(\lambda\) 是 \(A\) 的特征值，\(x\) 是 \(A\) 的特征向量。
那么 \(Ax=\lambda I x\)。\((\lambda I-A)x=0\)。因为 \(x\) 不为零向量，所以 \(\det(\lambda I-A)=0\)，也就是 \(\lambda I-A\) 不满秩。
我们称 \(\det(\lambda I-A)=0\) 为 \(A\) 的特征多项式，元是 \(\lambda\)。特征多项式是 \(n\) 次的，他的 \(n\) 个根就是 \(A\) 的所有特征值。（可能有相等的根）
对于上三角矩阵，所有的特征值就是主对角线上的所有值。
如果有 \(n\) 个线性无关的特征向量 \(x_i\)（当且仅当 \(A\) 满秩），那么有 \(A\begin{bmatrix}x_1&x_2&\cdots&x_n\end{bmatrix}=\begin{bmatrix}x_1&x_2&\cdots&x_n\end{bmatrix}\begin{bmatrix}\lambda_1&0&0&\cdots&0\\0&\lambda_2&0&\cdots&0\\\cdots\\0&0&0&\cdots&\lambda_n\end{bmatrix}\)。
Hamilton-Cayley 定理：对于矩阵 \(A\) 的特征多项式 \(f(\lambda)=\sum\limits_{i=0}^nc_i\lambda^i\)，有 \(f(A)=0\)，即 \(\sum\limits_{i=0}^nc_iA^i=0\)。
会了这些，就可以开始了。
\(O(k^3\log n)\) 的相信大家都会。（什么？不会？赶快去学矩阵快速幂）
首先考虑写出转移矩阵 \(A\) 和初始行向量 \(f\)，我们要求的是 \(f\times A^n\) 的第 \(0\) 维。
假如我们构造出了一个序列 \(c\) 使得
那么有：
那么就可以 \(O(k)\) 计算了。
那么 \(c\) 怎么弄呢？
令 \(R(A)=\sum\limits_{i=0}^{k-1}c_iA^i\)
假如存在 \(k\) 次多项式 \(G(A)\)，使得 \(A^n=F(A)G(A)+R(A)\)。(标准多项式除法形式)
当 \(G(A)=0\) 时就有 \(A^n=R(A)\)，所以要求的就是 \(A^n\bmod G(A)\)，快速幂+多项式除法 \(O(k\log k\log n)\) 解决。
那么如何构造 \(G(A)=0\) 的多项式呢？
看到上面的 Hamilton-Cayley 定理，令 \(G(\lambda)=\det(\lambda I-A)\) 即可。
手玩一下，发现 \(\det(\lambda I-A)=-\sum\limits_{i=0}^{k-1}\lambda^ia_{k-i}+\lambda^k\)。（注意交换行的时候取相反数！！！）
所以 \(g_i=-a_{k-i}(i\ne k),g_k=1\)。
时间复杂度 \(O(k\log k\log n)\)。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=333333,mod=998244353;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline int read(){
int x=0,f=0;char ch=getchar();
while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
return f?-x:x;
}
int n,k,s,f[maxn],g[maxn],fac[maxn],ans[maxn],prod[maxn],tmp[maxn],lim,l,rev[maxn],invt[maxn],Arev[maxn],Btmp[maxn],Brev[maxn],Brevinv[maxn],C[maxn];
inline int add(int x,int y){return x+y<mod?x+y:x+y-mod;}
inline int sub(int x,int y){return x<y?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int qpow(int a,int b){
int ans=1;
for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a);
return ans;
}
void init(int upr){
for(lim=1,l=0;lim<upr;lim<<=1,l++);
FOR(i,0,lim-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
}
void NTT(int *A,int tp){
FOR(i,0,lim-1) if(i<rev[i]) swap(A[i],A[rev[i]]);
for(int i=1;i<lim;i<<=1)
for(int j=0,r=i<<1,Wn=qpow(3,mod-1+tp*(mod-1)/r);j<lim;j+=r)
for(int k=0,w=1;k<i;k++,w=mul(w,Wn)){
int x=A[j+k],y=mul(w,A[i+j+k]);
A[j+k]=add(x,y);
A[i+j+k]=sub(x,y);
}
if(tp==-1){
int linv=qpow(lim,mod-2);
FOR(i,0,lim-1) A[i]=mul(A[i],linv);
}
}
void inv(int *A,int *B,int deg){
if(deg==1) return void(B[0]=qpow(A[0],mod-2));
inv(A,B,(deg+1)>>1);
init(deg<<1);
FOR(i,0,deg-1) invt[i]=A[i];
FOR(i,deg,lim-1) invt[i]=0;
NTT(invt,1);NTT(B,1);
FOR(i,0,lim-1) B[i]=mul(B[i],sub(2,mul(invt[i],B[i])));
NTT(B,-1);
FOR(i,deg,lim-1) B[i]=0;
}
void division(int *A,int *B,int *D,int n,int m){
if(n<m){
init(n<<1);
FOR(i,n+1,lim-1) D[i]=A[i]=0;
FOR(i,0,n) D[i]=A[i];
return;
}
init(n<<1);
FOR(i,0,lim-1) Arev[i]=Brev[i]=Brevinv[i]=C[i]=Btmp[i]=D[i]=0;
FOR(i,n+1,lim-1) A[i]=0;
FOR(i,m+1,lim-1) B[i]=0;
FOR(i,0,n) Arev[i]=A[n-i];
FOR(i,0,n-m) Brev[i]=B[m-i];
inv(Brev,Brevinv,n-m+1);
init(n<<1);
NTT(Arev,1);NTT(Brevinv,1);
FOR(i,0,lim-1) C[i]=mul(Arev[i],Brevinv[i]);
NTT(C,-1);
FOR(i,0,(n-m)/2) swap(C[i],C[n-m-i]);
FOR(i,n-m+1,lim-1) C[i]=0;
FOR(i,0,m) Btmp[i]=B[i];
NTT(Btmp,1);NTT(C,1);
FOR(i,0,lim-1) D[i]=mul(Btmp[i],C[i]);
NTT(D,-1);
FOR(i,0,m-1) D[i]=sub(A[i],D[i]);
FOR(i,m,lim-1) D[i]=0;
}
int main(){
n=read();k=read();
FOR(i,1,k) g[k-i]=(mod-read())%mod;
FOR(i,0,k-1) f[i]=(read()+mod)%mod;
g[k]=fac[1]=ans[0]=1;
while(n){
if(n&1){
init(k<<1);
NTT(fac,1);NTT(ans,1);
FOR(i,0,lim-1) prod[i]=mul(fac[i],ans[i]);
NTT(fac,-1);NTT(ans,-1);NTT(prod,-1);
division(prod,g,ans,(k-1)<<1,k);
}
init(k<<1);
NTT(fac,1);
FOR(i,0,lim-1) prod[i]=mul(fac[i],fac[i]);
NTT(fac,-1);NTT(prod,-1);
division(prod,g,fac,(k-1)<<1,k);
n>>=1;
}
FOR(i,0,k-1) s=add(s,mul(f[i],ans[i]));
printf("%d\n",s);
}
