常系数齐次线性递推
模板题：传送门
分析
严格规定 \(f_k\neq 0\)
不难列出转移矩阵：\(\left( \begin{matrix} a_{n} \\ a_{n-1} \\ a_{n-2} \\ \vdots \\ a_{n-k} \end{matrix} \right)= \left( \begin{matrix} f_1 & f_2 & \cdots & f_{k-1} & f_k \\ 1 \\ & 1 \\ && \ddots \\ &&& 1&0 \end{matrix} \right)\cdot \left( \begin{matrix} a_{n-1} \\ a_{n-2} \\ a_{n-3} \\ \vdots \\ a_{n-k-1} \end{matrix} \right)\)
记该转移为 \(\vec A_n=F_k\cdot \vec A_{n-1}\)
不难得出 \(\vec A_n=F_k^n \cdot \vec A_0\)
这里形式化定义了 $a_{-1}, a_{-2}, a_{-3},\cdots $
现在来进行一个求解：
假设我们能找到一个次数尽可能小的多项式函数 \(G\) 使得 \(G(F_k)=O_k\) 即 \(k\) 阶零矩阵
则我们可以得到： \(\vec A_n=F_k^n\cdot \vec A_0=[ Q(F_k)G(F_k)+R(F_k) ]\cdot \vec A_0=R(F_k)\cdot \vec A_0\)
其中，\(R\) 为 \(x^n \bmod G(x)\)
例如 \(n=2,G(x)=x^2-x-1\) 时，\(R(x)=x^2\bmod (x^2-x-1)=x+1\)
那么，显然有 \(\displaystyle \vec A_n=(\sum_{i=0}^{\text{deg }G-1} r_i\cdot F_k^i)\cdot \vec A_0=\sum_{i=0}^{\text{deg }G-1} r_i\cdot (F_k^i\cdot \vec A_0)=\sum_{i=0}^{\text{deg }G-1}r_i\cdot \vec A_i\)
\(\text{deg }G\) 表示多项式 \(G(x)\) 的次数
由于我们最后所求为 \(a_n\) ，仅为 \(\vec A_n\) 的第一个元素；故对等式右边向量 \(\vec A_i\) ，均只需取第一个元素 \(a_i\)
如果觉得不严谨，可以认为是方程两边同时右乘一个向量 \(\vec v=(1, 0, 0, \cdots , 0)\)
故 \(\displaystyle a_n=\sum_{i=0}^{\text{deg }G-1}r_i\cdot a_i\)
现在我们考虑如何构造 \(G(x)\) 使得 \(G(F_k)=O_k\)
先给出结论，后续提供证明：
\(\text{deg }G=k, g_n=\begin{cases} -f_{k-n}, n<k \\ 1, n=k \end{cases}\)
故原式变换为 \(\displaystyle a_n=\sum_{i=0}^{k-1} r_i\cdot a_i\)
其中，\(a_0\)~\(a_{k-1}\) 为题目给定的初值，\(r_0\)~\(r_{k-1}\) 为 \(x^n\bmod G(x)\) 的多项式取余
对于求解 \(x^n\bmod G(x)\) 可以由倍增实现，故复杂度为
无 FFT 优化：\(O(k^2\log n)\)；有 FFT 优化：\(O(k\log k\cdot \log n)\)
附 1 ：AC 代码
FFT 版：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
typedef pair<int,int> pii;
#define fi first
#define se second
#define de(x) cout<<#x<<" = "<<x<<endl
#define dd(x) cout<<#x<<" = "<<x<<" "
const int LimBit=18, M=1<<LimBit<<1, MOD=998244353;
int a[M], b[M], c[M];
inline int Normal(int n) { return (n%=MOD)>=0?n:n+MOD; }
struct NTT{
int G=3, P=998244353;
int N, na, nb, W[2][M+M], Rev[M+M], *w[2], *rev, invN, InV[M];
inline ll kpow(ll a, ll x) { ll ans=1; for(;x;x>>=1,a=a*a%P) if(x&1) ans=ans*a%P; return ans; }
inline ll inv(ll a) { return kpow(a, P-2); }
inline int add(int a, int b) { return (a+b>=P)?(a+b-P):(a+b); }
inline int dis(int a, int b) { return (a-b<0)?(a-b+P):(a-b); }
inline void display(int *f,int n){
for(int i=0;i<n;++i) cout<<f[i]<<" "; cout<<"\n";
}
NTT(){
InV[1] = 1;
for(int i=2;i<P&&i<M;++i)
InV[i]=P-(ll)P/i*InV[P%i]%P;
w[0]=W[0], w[1]=W[1], rev=Rev;
w[0][0]=w[1][0]=1;
for(int i=1, x=kpow(G,P>>LimBit), y=inv(x); !(i>>LimBit); ++i){
rev[i]=(rev[i>>1]>>1)|((i&1)<<LimBit-1);
w[0][i]=(ll)w[0][i-1]*x%P, w[1][i]=(ll)w[1][i-1]*y%P;
}
int *lstw[2]={w[0], w[1]};
w[0]+=1<<LimBit, w[1]+=1<<LimBit, rev+=1<<LimBit;
for(int d=LimBit-1; d>=0; --d){
for(int i=0; !(i>>d); ++i){
rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
w[0][i]=lstw[0][i<<1], w[1][i]=lstw[1][i<<1];
}
lstw[0]=w[0], lstw[1]=w[1];
w[0]+=1<<d, w[1]+=1<<d, rev+=1<<d;
}
}
inline void work(){
w[0]=W[0], w[1]=W[1], rev=Rev;
for(int d=LimBit;1<<d>N;--d)
w[0]+=1<<d, w[1]+=1<<d, rev+=1<<d;
invN=inv(N);
}
inline void FFT(int *a,int f){
for(int i=0;i<N;++i) if(i<rev[i]) swap(a[i], a[rev[i]]);
for(int i=1;i<N;i<<=1)
for(int j=0, t=N/(i<<1); j<N; j+=i<<1)
for(int k=0, l=0, x; k<i; ++k, l+=t)
x=(ll)w[f][l]*a[j+k+i]%P, a[j+k+i]=dis(a[j+k], x), a[j+k]=add(a[j+k], x);
if(f) for(int i=0;i<N;++i) a[i]=(ll)a[i]*invN%P;
}
inline void doit(int *a,int *b,int na,int nb,int *c=0){
static int x[M], y[M];
if(!c) c=a;
for(N=1;N<na+nb-1;N<<=1);
memset(x, 0, sizeof(x[0])*N); memset(y, 0, sizeof(y[0])*N);
memcpy(x, a, sizeof(a[0])*na); memcpy(y, b, sizeof(b[0])*nb);
work(); FFT(x, 0); FFT(y, 0);
for(int i=0;i<N;++i) x[i]=(ll)x[i]*y[i]%P;
FFT(x, 1);
memcpy(c, x, sizeof(x[0])*N);
}
inline void get_inv(int *f,int *g,int n){
g[0]=inv(f[0]);
if(n==1) return ;
for(int l=0; (1<<l)<n; ++l){
memcpy(a, f, sizeof(a[0])<<l+1);
memset(a+(1<<l+1), 0, sizeof(a[0])<<l+1);
memset(g+(1<<l), 0, sizeof(g[0])*3<<l);
N=1<<l+2;
work(); FFT(a, 0); FFT(g, 0);
for(int i=0; i<N; ++i) g[i]=(ll)g[i]*dis(2, (ll)g[i]*a[i]%P)%P;
FFT(g, 1);
}
memset(g+n, 0, sizeof(g[0])*(N-n));
}
inline void get_div(int *f,int *g,int *q,int *r, int n,int m){
while(f[n-1]==0&&n) --n;
while(g[m-1]==0&&m) --m;
if(n<m||m==0){
memcpy(r, f, sizeof(f[0])*n);
memset(q, 0, sizeof(q[0])*n);
return ;
}
reverse(g, g+m);
get_inv(g, q, n-m+1);
reverse(f, f+n);
doit(q, f, n-m+1, n-m+1);
reverse(q, q+n-m+1);
reverse(g, g+m);
reverse(f, f+n);
memset(q+n-m+1, 0, sizeof(q[0])*(m-1) );
memcpy(r, g, sizeof(g[0])*m);
doit(r, q, m, n);
for(int i=0;i<m-1;++i) r[i]=dis(f[i], r[i]);
}
inline void get_xpow(int *g,int *m,int x,int nm){
if(x==0){
g[0]=1;
return ;
}
get_xpow(g, m, x>>1, nm);
for(N=1;N<nm+nm-3;N<<=1);
for(int i=nm-1;i<N;++i) g[i]=0;
work(); FFT(g, 0);
for(int i=0;i<N;++i) g[i]=(ll)g[i]*g[i]%MOD;
FFT(g, 1);
if(x&1){
for(int i=N;i>=1;--i) g[i]=g[i-1];
g[0]=0;
get_div(g, m, c, b, N+1, nm);
}
else get_div(g, m, c, b, N, nm);
for(N=1;N<nm+nm-3;N<<=1);
for(int i=nm-1;i<N;++i) g[i]=0;
for(int i=0;i<nm-1;++i) g[i]=b[i];
}
}ntt;
int n, k, f[M], g[M], r[M];
inline int ans(){
int res=0;
for(int i=0;i<k;++i) res=ntt.add(res, (ll)g[i]*r[i]%MOD);
return res;
}
inline void init(){
cin>>n>>k;
for(int i=1;i<=k;++i) cin>>f[k-i], f[k-i]=Normal(-f[k-i]);
for(int i=0;i<k;++i) cin>>g[i], g[i]=Normal(g[i]);
f[k]=1;
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
init();
ntt.get_xpow(r, f, n, k+1);
cout<<ans();
// cout<<fixed<<setprecision(6);
cout.flush();
return 0;
}
附 2 ：多项式取余的理论
设给定 \(F(x), G(x)\) ，求解 \(F(x)=Q(x)\cdot G(x)+R(X)\)
其中，\(Q(x)\) 为多项式 \(F(x)\) 除以 \(G(x)\) 的商，\(R(x)\) 是除法的余数
观察可以发现一个很棒的性质：当 \(\text{deg }F\geq \text{deg }G\) 时，有 \(\text{deg }Q=\text{deg }F-\text{deg G},\text{deg }R= \text{deg }G-1\)
这里的 \(R(x)\) 强制设置为最高可能的位数
故我们考虑将 \(x\) 换元为 \({1\over x}\) 且在方程两边同时乘上 \(x^{\text{deg } F}\) 得到：
\(x^{\text{deg }F}F({1\over x})=x^{\text{deg }F-\text{deg }G}Q({1\over x})\cdot x^{\text{deg }G}G({1\over x})+x^{\text{deg }F-G+1}\cdot x^{\text{deg }G-1}R({1\over x})\)
让我们再考虑一下 \(x^{\text{deg }F}F({1\over x})\) 的含义：
原来的 \(n\) 次项 \(f_n\cdot x^n\) 转变为了 \(f_n\cdot {1\over x^n}\cdot x^{\text{deg }F}=f_n\cdot x^{\text{deg }F-n}\) 成为了现在的 \((\text{deg }F-n)\) 次项
相当于原来的多项式函数 \(F(x)\) 的后 \(\text{deg }F\) 次项镜像翻转了一下，故我们记这种多项式为 \(F^R(x)\)
代入上面得到 \(F^R(x)=Q^R(x)\cdot G^R(x)+x^{\text{deg }F-\text{deg }G+1}\cdot R^R(x)\)
两边同余 \(x^{\text{deg }F-\text{deg }G+1}\) 得 \(F^R(x)\equiv Q^R(x)\cdot G^R(x)\pmod{x^{\text{deg }F-\text{deg }G+1}}\)
移项得到 \(Q^R(x)\equiv [F^R(x)]^{-1}\cdot G^R(x)\pmod{x^{\text{deg }F-\text{deg }G+1}}\)
通过多项式的取逆，我们可以得到 \(Q^R(x)\) 的后 \((\text{deg }F-\text{deg }G)\) 项；不过因为 \(\text{deg }Q=\text{deg }F-\text{deg }G\) ，我们直接得到了 \(Q\) 的所有项
故两多项式乘积后，直接翻转得到 \(Q(x)\)
再通过多项式乘法和多项式减法求得 \(R(x)=F(x)-Q(x)\cdot G(x)\)，复杂度为 \(O(n\log n)\)
附 3 ： \(G(x)\) 的构造
根据 Cayley-Hamilton 定理，对于 \(F_k\) 的特征多项式 \(p(\lambda)\) 有 \(p(F_k)=O_k\) 即 \(k\) 阶零矩阵
故只要令 \(G(x)\) 为 \(F_k\) 的特征多项式即可
现考虑如何求解 \(F_k\) 的特征多项式：
\(p(\lambda)=\det (F_k-\lambda\cdot E_k)=\det \left( \begin{matrix} f_1-\lambda & f_2 & \cdots & f_{k-1} & f_k \\ 1 & -\lambda \\ & 1 & -\lambda \\ && \ddots & \ddots \\ &&& 1&-\lambda \end{matrix} \right)\)
据说有人手动展开后，得到 \(p(\lambda)=(-1)^n\cdot (\lambda^n-f_1\cdot \lambda^{n-1}-f_2\cdot \lambda^{n-2}-\cdots -f_k\cdot \lambda^0)\)
由于 \(G(F_k)=O_k\) ，所以前面那个 \((-1)^n\) 也不怎么重要了，直接得到：
\(\displaystyle G(x)=x^n-\sum_{i=1}^kf_ix^{n-i}\)
于是就有了上面那个式子
拓展
有个神仙算法叫 Berlekamp-Massey 算法
将某个序列的前若干项丢进去，它会贪心地得出满足这些项的最短递推式
设丢进去的项数为 \(k\) ，则复杂度为 \(O(k^2)\)
结合该算法，可以得出一个非常优（暴）美（力）的算法：
先暴力打表，尽可能多打几项，然后丢进 BM 算法，跑出最短递推式，然后用常系数齐次线性递推算法直接得出第 \(n\) 项
由于最短递推式最多为 \(k\) 项，故总复杂度为 \(O(k^2+k\log k\log n)\) 或 \(O(k^2\log n)\)
相比于拉格朗日插值，优势非常明显：
对于类似卡特兰数 \(\displaystyle c_n={1\over n+1}\dbinom{2n}{n}\) ，其线性项数不定，拉格朗日插值无从下手
但其最短递推式并不长：\(\displaystyle c_n={1\over n+1}\dbinom{2n}{n}={(2n)!\over n!(n+1)!}={(2n)\cdot (2n-1)\over (n+1)\cdot n}\cdot {(2n-2)!\over (n-1)!n!}={4n-2\over n+1}\cdot {1\over n}\dbinom{2n-2}{n-1}={4n-2\over n+1}c_{n-1}\)
一项的递推式，显然跑得飞快
板子的话，显然网上有，我就不弄了
题意：给定一条 \(n\) 元链，链上相邻元素相连，共有 \((n-1)\) 条边。链外有 \(3\) 个点，每个点连向链上的 \(n\) 个点。整个图共 \((4n-1)\) 条边。\(n\leq 10^9\) ，问生成树方案数，答案对 \((10^9+7)\) 取模。
当图中点数为 \(|V|\) 时，根据矩阵树定理：可以建立基尔霍夫矩阵，用高斯消元法 \(O(|V|^3)\) 求出生成树的方案数。
由于点数为 \(n+3\) ，故我们可以暴力打出约前 \(100\) 项的生成树方案数
接着，丢进 BM 算法，跑出最短递推式，发现只有 \(6\) 项，为：
\(\displaystyle a_n=\sum_{i=1}^6 f_ia_{n-i}\)
其中，\(f=\{-1, 15, -78, 155, -78, 15\}\)
直接用常系数齐次线性递推求出答案即可
