【学习笔记】NTT 与多项式进阶操作

NTT

在模 (p) 意义下，若 (p) 的原根为 (g)，可以用 (g^{(p-1)/n}) 代替 (omega_n^1)，它满足 FFT 中单位根的所有性质。

但必须 (nmid p-1) 才行。我们通常把 (n) 补为 (2) 的幂。所以如果 (p=r imes 2^m+1)，可以处理项数为 (2^m) 的多项式。

const int N = (1 << 20) + 1, P = 998244353, G = 3;
ll invG;
int tr[N], n;

inline ll poww(ll x, ll y = P - 2) {
	ll res = 1;
	for(; y; y >>= 1, x = x * x % P)
		if(y & 1) res = res * x % P;
	return res;
}

inline void madd(ll &x, const ll &y) {
	x += y;
	if(x >= P) x -= P;
}

inline ll mod(const ll x) {
	return x >= P ? x - P : x;
}

inline void cpy(ll *f, ll *g, int n) {
	rep(i, 0, n - 1) *(f + i) = *(g + i);
}

inline void clear(ll *f, int n) {
	rep(i, 0, n - 1) *(f + i) = 0;
}

inline void mul(ll *f, ll *g, int n) {
	rep(i, 0, n - 1) (f[i] *= g[i]) %= P;
}

void NTT(ll *f, int n, int fl) { // fl=1表示DFT，fl=0表示IDFT
	ll tmp, w1, buf;
	rep(i, 0, n - 1) tr[i] = (tr[i >> 1] >> 1) | (i & 1 ? (n >> 1) : 0);
	rep(i, 0, n - 1) if(i < tr[i]) swap(f[i], f[tr[i]]);
	for(rei lim = 2, lim2; lim <= n; lim <<= 1) {
		w1 = poww(fl ? G : invG, (P - 1) / lim);
		lim2 = lim >> 1;
		for(rei i = 0; i < n; i += lim) {
			buf = 1;
			rep(k, i, i + lim2 - 1) {
				tmp = buf * f[k + lim2] % P;
				f[k + lim2] = mod(f[k] - tmp + P);
				madd(f[k], tmp);
				(buf *= w1) %= P;
			}
		}
	}
	if(!fl) {
		tmp = poww(n);
		rep(i, 0, n) (f[i] *= tmp) %= P;
	}
}

1. 乘法逆元

题目链接

若 (F(x)) 为 (n-1) 次多项式，(F(x)G(x)equiv 1~(mod~x^n))，则称 (G(x)) 是 (F(x)) 的乘法逆元。

如果 (F(x)) 有乘法逆元，则 (F[0] eq 0) （(F[i]) 表示 (F(x)) 的 (i) 次项系数）。因为如果 (F[0]=0)，无论 (G[0]) 是否为 (0)，(F(x)G(x)) 的常数项都为 (0)。

在模 (x) 意义下，(F(x)) 的乘法逆元为 (F[0]^{-1})。

考虑我们已经知道 (F(x)) 在 (mod~x^{n/2}) 意义下，逆元为 (G^prime(x)) 即 (F(x)G^prime(x)equiv 1~(mod~x^{n/2}))。

现求满足 (F(x)G(x)equiv 1~(mod~x^n)) 的 (G(x))。两式相减可得 (F(x)(G(x)-G^prime(x))equiv 0~(mod~x^{n/2}))

显然 (F(x) otequiv 0~(mod~x^{n/2}))，那么 ((G(x)-G^prime(x))equiv 0~(mod~x^{n/2}))。平方可得 (G(x)^2-2G(x)G^prime(x)+G^prime(x)^2equiv 0~(mod~x^n))。

两边同乘 (F(x)) 得 (G(x)^2F(x)-2F(x)G(x)G^prime(x)+F(x)G^prime(x)^2equiv 0~(mod~x^n)) 即 (G(x)-2G^prime(x)+F(x)G^prime(x)^2equiv 0~(mod~x^n))。

所以 (G(x)equiv 2G^prime(x)-F(x)G^prime(x)^2~(mod~x^n))。

倍增求解即可。

复杂度：(T(n)=T(dfrac{n}{2})+operatorname{O}(nlog n)=operatorname{O}(nlog n))。

void poly_inv(ll *f, int m) {
	static int n;
	static ll s[N], h[N], g[N];
	for(n = 1; n < m; n <<= 1) ;
	clear(s, n * 2);
	clear(g, n * 2);
	clear(h, n * 2);
	g[0] = poww(f[0]);
	for(rei lim = 2; lim <= n; lim <<= 1) {
		rep(i, 0, lim / 2 - 1) h[i] = mod(g[i] * 2);
		cpy(s, f, lim);
		NTT(g, lim * 2, 1);
		mul(g, g, lim * 2);
		NTT(s, lim * 2, 1);
		mul(g, s, lim * 2);
		NTT(g, lim * 2, 0);
		rep(i, 0, lim - 1) g[i] = mod(h[i] - g[i] + P);
		clear(g + lim, lim); // 一定要清空！！！！！！
	}
	cpy(f, g, m);
}

2. 开根

题目链接

在模 (x) 意义下，直接求二次剩余。本题保证 (a_0=1)，就不需要求了，但加强版需要。

考虑我们已经知道在 (mod~x^{n/2}) 意义下，有 (G^prime(x)^2equiv F(x)~(mod~x^{n/2}))。

现求满足 (G(x)^2equiv F(x)~(mod~x^n)) 的 (G(x))。

显然有 (G(x)equiv G^prime(x)~(mod~x^{n/2}))。

移项，平方，得 (G(x)^2-2G(x)G^prime(x)+G^prime(x)^2equiv 0~(mod~x^n)) 即 (F(x)-2G(x)G^prime(x)+G^prime(x)^2equiv 0~(mod~x^n))

整理得 (G(x)equiv dfrac{1}{2}(G^prime(x)+dfrac{F(x)}{G^prime(x)}))。

倍增求解+求逆即可。

复杂度：(T(n)=T(dfrac{n}{2})+operatorname{O}(nlog n)=operatorname{O}(nlog n))。

void poly_sqrt(ll *f, int m) {
	static ll g[N], h[N], s[N];
	static int n;
	for(n = 1; n < m; n <<= 1) ;
	clear(s, n * 2);
	clear(g, n * 2);
	clear(h, n * 2);
	g[0] = 1; //由于保证a[0]=1，1^2=1。如果没有保证需求二次剩余
	for(rei lim = 2; lim <= n; lim <<= 1) {
		rep(i, 0, lim / 2 - 1) h[i] = g[i];
		clear(h + lim / 2, lim / 2);
		poly_inv(h, lim);
		NTT(h, lim * 2, 1);
		cpy(s, f, lim);
		NTT(s, lim * 2, 1);
		mul(h, s, lim * 2);
		NTT(h, lim * 2, 0);
		rep(i, 0, lim - 1) g[i] = (g[i] + h[i]) * inv2 % P;
		clear(g + lim, lim);
	}
	cpy(f, g, m);
}

3. 对数函数（ln)

题目链接

导数

百度百科

函数在某一点的导数就是该函数所代表的曲线在这一点上的切线斜率。

公式：

(y=a,y^prime=0) (C为常数）

(y=x^n,y^prime=nx^{n-1})

(y=ln x,y^prime=dfrac{1}{x})

(f(g(x))^prime=f^prime(g(x))g^prime(x))

不定积分

百度百科

求导的逆运算。

(int a;dx=ax) (a为常数）

(int x^n;dx=dfrac{x^{n+1}}{n+1})

(int dfrac{1}{x};dx=ln vert xvert)

本题

下文中省略 ((mod~998244353))。

(G(x)equivln(F(x)))

(G^prime(x)equiv ln(F(x))^primeequiv ln^prime(F(x))F^prime(x)equivdfrac{F^prime(x)}{F(x)})

(G(x)equivintdfrac{F^prime(x)}{F(x)};dx)

void poly_der(ll *f, int n) {
	rep(i, 1, n - 1) f[i - 1] = f[i] * i % P;
	f[n - 1] = 0;
}

void poly_int(ll *f, int n) {
	per(i, n - 1, 0) f[i + 1] = f[i] * inv[i + 1] % P;
	f[0] = 0;
}

void poly_ln(ll *f, int m) {
	static ll g[N];
	static int n;
	for(n = 1; n < m * 2; n <<= 1) ;
	cpy(g, f, m);
	poly_inv(g, m);
	poly_der(f, m);
	NTT(g, n, 1);
	NTT(f, n, 1);
	mul(f, g, n);
	NTT(f, n, 0);
	poly_int(f, m);
}

4. 指数函数（exp）

题目链接

百度百科 - 牛顿迭代

考虑使用牛顿迭代求一个函数的零点，若某一次得到结果 (x)，则下一次结果为 (x-dfrac{f(x)}{f^prime(x)})

多项式同理。经验可得 每迭代一次，精度会从 (x^k) 变成 (x^{2k})。

即：若 (F(G_0(x))=0~(mod~x^{n/2})) 且 (F(G(x))=0~(mod~x^n))，则 (G(x)=G_0(x)-dfrac{F(G_0(x))}{F^prime(G_0(x))})。

证明

令 (G(x)=e^{F(x)})

(ln G(x)-F(x)=0)

如果等号左边是一个关于 (G(x)) 的函数 (H(G(x)))，把 (F(x)) 当作常数项，就可以牛顿迭代做。

(G(x)=G_0(x)-dfrac{H(G_0(x))}{H(G_0(x))^prime}=G_0(x)-dfrac{ln G_0(x)-F(x)}{dfrac{1}{G_0(x)}}=G_0(x)(1-ln G_0(x)+F(x)))

倍增求解即可。

void poly_exp(ll *f, int m) {
	static ll g[N], h[N];
	static int n;
	for(n = 1; n < m; n <<= 1) ;
	clear(g, n * 2); clear(h, n * 2);
	g[0] = 1;
	for(rei lim = 2; lim <= n; lim <<= 1) {
		cpy(h, g, lim / 2);
		poly_ln(h, lim);
		rep(i, 0, lim - 1) h[i] = mod(f[i] - h[i] + P);
		madd(h[0], 1);
		NTT(g, lim * 2, 1);
		NTT(h, lim * 2, 1);
		mul(g, h, lim * 2);
		NTT(g, lim * 2, 0);
		clear(g + lim, lim);
		clear(h + lim, lim);
	}
	cpy(f, g, m);
}

5.快速幂

题目链接

若保证 (f[0]=1)，可以

(G(x)equiv F^k(x))

(ln G(x)equiv kln F(x))

(G(x)equiv e^{kln F(x)})

若 (f[0]>1)，可以将整个多项式除以 (f[0])

若 (f[0]=0)，找到系数不为 (0) 的最低次项，进行一系列位移即可

6. 带余除法

题目链接

可以推导 (F(x)) 与 (F(frac{1}{x})) 的关系：(F(frac{1}{x})=sumlimits_{i=0}^nf_i(frac{1}{x})^i=dfrac{sumlimits_{i=0}^nf_ix^{n-i}}{x^n}=dfrac{F^*(x)}{x^n})（(F^*(x)) 是将 (F(x)) 倒过来）

(F(x)=Q(x)G(x)+R(x))

(F(frac{1}{x})=Q(frac{1}{x})G(frac{1}{x})+R(frac{1}{x}))

(F^*(x)=Q^*(x)G^*(x)+x^{n-m+1}R^*(x))

对 (x^{n-m+1}) 取模

(F^*(x)equiv Q^*(x)G^*(x)~(mod x^{n-m+1}))

可求得 (Q^*(x))（注意模数）

然后 (R(x)=F(x)-Q(x)G(x))