从傅里叶变换(FFT)到数论变换(NTT)

FFT可以用来计算多项式乘法，但是复数的运算中含有大量的浮点数，精度较低。对于只有整数参与运算的多项式，有时，( ext{NTT(Number-Theoretic Transform)})会是更好的选择。

阶

若(a,p)互素，且(p>1)，对于(a^k equiv 1 (mod p))的最小的(k)，称为(a)模(p)的阶，记做(sigma_p(a))。

(E.g.) (sigma_7(2)=3)

(2^1equiv 2(mod 7))

(2^2equiv 4(mod 7))

(2^3equiv 1(mod 7))

对于一个数(g)，(g)的阶一定是(p-1)的约数。

证明：

假设最小的(k)不是(p-1)的约数，找到(x)满足(xk>p-1>(x-1)k)，由费马小定理可知

[g^{xk}equiv g^{p-1}equiv 1 equiv g^{xk-(p-1)} (mod p) ]

(xk-(p-1)<k)，与假设矛盾。

原根

定义

(FFT)中，我们使用单位复根(omega_n^k=cos kfrac{2pi}{n}+isin kfrac{2pi}{n})，那有没有什么能够代替单位复根且解决精度问题呢？这就是原根。

设(m)是正整数，(a)是整数，若(a)模(m)的阶等于(varphi(m))，则称(a)为模(m)的一个原根。

若(p)为质数，设(g)为(p)的原根，那么(g^i mod p(1<j<p,1le ile p-1))的结果两两不同。且其等价于(g^{p-1}equiv 1(mod p))当且仅当指数为(p-1)的时候成立。(这里(p)是素数)

简单证明一下：

显然(g^0 equiv 1(mod p))

由原根的定义可知满足(g^{i} equiv 1(mod p))的最小正整数为(varphi(p)=p-1)

故由指数循环节可知，(g^i mod p(1<j<p,1le ile p-1))的结果两两不同。

性质

考虑在FFT当中我们需要单位复根的以下性质：

1. (omega_n^t)互不相同，保证点值的合法性；

2. (omega_{2n}^{2k} = omega_n^k)，用于分治；

3. (omega_n^{k+frac{n}{2}} = -omega_n^k)，用于分治；

4. 当(k eq 0)时，(1+omega_n^k+(omega_n^k)^2+dots +(omega_n^k)^{n-1}=0)，用于逆变换。

性质一

令(omega_n=g^q)，(1,g^q,g^{2q},dots,g^{(n-1)q})互不相同，满足性质一。

性质二

由(omega_n = g^q)可知，(omega_{2n}=g^{frac{q}{2}}(p=frac{q}{2} imes 2n + 1))，故(omega_{2n}^{2k} = g^{2k frac{q}{2}}=g^{kq}=g^q)，满足性质二。

性质三

根据费马小定理得

[omega_n^n=g^{nq}=g^{p-1}equiv 1(mod p) ]

又因为((omega_n^{frac{n}{2}})^2=omega_n^n)，所以(omega_n^{frac{n}{2}}equiv pm 1 (mod p))，而根据性质一可得(omega_n^{frac{n}{2}} eq omega_n^0)，即(omega_n^{frac{n}{2}}equiv -1(mod p))。可推出(omega_n^{k+frac{n}{2}}=omega_n^k imes omega_n^{frac{n}{2}}=-omega_n^k (mod p))，满足性质三。

性质四

当(k eq 0)时

[S(omega_n^k)=1+omega_n^k+(omega_n^k)^2+dots +(omega_n^k)^{n-1} ]

[omega_n^k S(omega_n^k)=omega_n^k+(omega_n^k)^2+(omega_n^k)^3+dots +(omega_n^k)^n ]

[omega_n^k S(omega_n^k)-S(omega_n^k)=(omega_n^k)^n-1 ]

[S(omega_n^k)=frac{(omega_n^k)^n-1}{omega_n^k-1} ]

由性质三中的推论可知，((omega_n^k)^n-1=(omega_n^n)^k-1equiv omega_n^n-1equiv 0 (mod p))，故(S(omega_n^k)=0)，性质四成立。

求原根

求一个质数的原根，可以用枚举法——枚举(g)，检验(g)是否为(p)的原根。

根据前面的关于阶知识可知，检验时，只需枚举(p-1)的所有约数(q)，检验(g^qequiv 1(mod p))即可。

代码实现

将(FFT)里所有关于(omega_n)的运算替换成(g^q)在模意义下的运算即可，注意(div n)要改为( imes n^{-1})

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
inline ll ty() {
	char ch = getchar(); ll x = 0, f = 1;
	while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
	while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
	return x * f;
}

const int _ = 4e6 + 10;
const ll P = 998244353, G = 3, Gx = 332748118;
int N, M, r[_];
ll A[_], B[_];

ll ksm(ll a, ll b) {
	ll ret = 1;
	for ( ; b; b >>= 1) {
		if (b & 1) ret = ret * a % P;
		a = a * a % P;
	}
	return ret;
}

void ntt(int lim, ll *a, int op) {
	for (int i = 0; i < lim; ++i)
		if (i < r[i]) swap(a[i], a[r[i]]);
	for (int len = 2; len <= lim; len <<= 1) {
		int mid = len >> 1;
		ll Wn = ksm(op == 1 ? G : Gx, (P - 1) / len);
		for (int i = 0; i < lim; i += len) {
			ll w = 1;
			for (int j = 0; j < mid; ++j, w = (w * Wn) % P) {
				ll x = a[i + j], y = w * a[i + j + mid] % P;
				a[i + j] = (x + y) % P;
				a[i + j + mid] = (x - y + P) % P;
			}
		}
	}
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("ntt.in", "r", stdin);
	freopen("ntt.out", "w", stdout);
#endif
	N = ty(), M = ty();
	for (int i = 0; i <= N; ++i) A[i] = (ty() + P) % P;
	for (int i = 0; i <= M; ++i) B[i] = (ty() + P) % P;
	int lim = 1, k = 0;
	while (lim <= N + M) lim <<= 1, ++k;
	for (int i = 0; i < lim ; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
	ntt(lim, A, 1);
	ntt(lim, B, 1);
	for (int i = 0; i < lim; ++i) A[i] = (A[i] * B[i]) % P;
	ntt(lim, A, -1);
	ll invx = ksm(lim, P - 2);
	for (int i = 0; i <= N + M; ++i)
		printf("%lld ", (A[i] * invx) % P);
	return 0;
}

参考资料

从傅里叶变换到数论变换 | Menci's Blog

快速数论变换(NTT)小结 - 自为风月马前卒 - 博客园