一些閑話 & 免責聲明：感覺好久沒有更過學習筆記了，現在才來學這個不知道來不來得及呢。因為筆者是一個數學菜雞，所以并不會很嚴謹（。
問題描述
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;
}
浙公網安備 33010602011771號