[学习笔记] 常系数齐次线性递推

一些闲话 & 免责声明：感觉好久没有更过学习笔记了，现在才来学这个不知道来不来得及呢。因为笔者是一个数学菜鸡，所以并不会很严谨（。

问题描述

\(\scr Portal.\)

LSB-first 算法

从线性递推式到分式第 \(n\) 项

记 \(Q(x)=1-f_1x-f_2x^2-\dots -f_kx^k\)，\(F ( x )\) 为递推数列的 \(\tt OGF\)。我们需要求解 \(F(x)\) 的 \(n\) 次项系数。这里直接给出一个结论：

存在一个最高次项的次数小于 \(k\) 的多项式 \(P(x)\) 满足

\[F(x)=\frac{P(x)}{Q(x)} \]

由于 \(P(x)\) 的度数 \(<k\)，而 \(F(x)\) 的前 \(k\) 项是已知的，所以可以仅用 \(F(x)\) 的前 \(k\) 项与 \(Q ( x )\) 进行多项式乘法求出 \(P ( x )\).

剩下的问题是，要怎么快速计算 \(\displaystyle \frac{P(x)}{Q(x)}\) 的第 \(n\) 项。

Bostan-mori 算法

首先有

\[[x^n]\frac{P(x)}{Q(x)}=[x^n]\frac{P(x)Q(-x)}{Q(x)Q(-x)} \]

可以发现，\(Q(x)Q(-x)\) 是一个偶函数（考虑将 \(x,-x\) 分别代入），也就是说，它的奇次项均为零，于是可以构造多项式 \(v(x^2)=Q(x)Q(-x)\)，同时对分子奇偶分类可得

\[\begin{align}&[x^n]\frac{c_{\rm even}(x^2)+x\cdot c_{\rm odd}(x^2)}{v(x^2)}\\=&\begin{cases} [x^{n/2}]\frac{c_{\rm even}(x)}{v(x)}, & \text{if }n\text{ is even}\\ [x^{n/2}]\frac{c_{\rm odd}(x)}{v(x)}, & \text{if }n\text{ is odd}\end{cases}\end{align} \]

对于后面的一步推导，考虑 \(v(x)\) 的逆也是偶函数即可。于是可以进行递推求解。复杂度 \(\mathcal O(k\log k\log n)\).

代码实现

三个 \(998244353\) 相加又忘记开 long long 了，我自裁 🔫。

/* 先复习动人的墓志，再背诵浪漫的情诗。 */

# include <cstdio>
# include <cctype>
# define print(x,y) write(x), putchar(y)

template <class T>
inline T read(const T sample) {
    T x=0; char s; bool f=0;
    while(!isdigit(s=getchar())) f|=(s=='-');
    for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
    return f? -x: x;
}
template <class T>
inline void write(T x) {
    static int writ[50], w_tp=0;
    if(x<0) putchar('-'), x=-x;
    do writ[++w_tp]=x-x/10*10, x/=10; while(x);
    while(putchar(writ[w_tp--]^48), w_tp);
}

# include <cstring>
# include <iostream>
using namespace std;

const int MAXN = 170000;
const int mod = 998244353, G = 3;

void mul(int& x,int y) { x=1ll*x*y%mod; }
int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
int inv(int x,int y=mod-2,int r=1) {
    for(; y; y>>=1, x=1ll*x*x%mod)
        if(y&1) r=1ll*r*x%mod; return r;
} const int iG = inv(G);

int P[MAXN], rev[MAXN], g[MAXN];
int lim, bit, Inv, n, k, Q[MAXN];

void ntt(int* f,bool opt=0) {
    for(int i=0;i<lim;++i)
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1) {
        const int wn = inv(opt?iG:G,(mod-1)/(mid<<1));
        for(int i=0; i<lim; i+=(mid<<1)) { int w=1;
            for(int j=0; j<mid; ++j, w=1ll*w*wn%mod) {
                const int tmp = 1ll*f[i|j|mid]*w%mod;
                f[i|j|mid] = dec(f[i|j], tmp),
                f[i|j] = inc(f[i|j], tmp);
            }
        }
    } if(!opt) return;
    for(int i=0;i<lim;++i) mul(f[i],Inv);
}

int main() {
    n=read(9), k=read(9); Q[0]=lim=1;
    for(int i=1;i<=k;++i) Q[i] = (0ll+(mod<<1)-read(9))%mod;
    for(int i=0;i<k;++i) P[i] = (0ll+(mod<<1)+read(9))%mod;
    while(lim<(k<<1|1)) lim<<=1, ++bit; 
    for(int i=0;i<lim;++i)
        rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1)); 
    ntt(P), ntt(Q); Inv = inv(lim);
    for(int i=0;i<lim;++i) mul(P[i],Q[i]);
    ntt(P,1), ntt(Q,1);
    const int SIZE = sizeof(int)*(lim-k-1);
    memset(P+k,0,sizeof(int)*(lim-k));
    memset(Q+k+1,0,SIZE);
    for(; n; n>>=1) {
        for(int i=0;i<=k;++i)
            g[i] = (i&1)?(mod-Q[i])%mod:Q[i];
        memset(g+k+1,0,SIZE);
        ntt(g), ntt(P), ntt(Q);
        for(int i=0;i<lim;++i)
            mul(P[i],g[i]), mul(Q[i],g[i]);
        ntt(P,1), ntt(Q,1);
        for(int i=0;i<=k;++i) Q[i]=Q[i<<1];
        for(int i=0;i<=k;++i) P[i]=P[i<<1|(n&1)];
        memset(P+k+1,0,SIZE), memset(Q+k+1,0,SIZE);
    } print((mul(P[0],inv(Q[0])), P[0]), '\n');
    return 0;                               
}