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

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

FFT可以用来计算多项式乘法，但是复数的运算中含有大量的浮点数，精度较低。对于只有整数参与运算的多项式，有时，\(\text{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^1\equiv 2(\mod 7)\)

\(2^2\equiv 4(\mod 7)\)

\(2^3\equiv 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 k\frac{2\pi}{n}+i\sin k\frac{2\pi}{n}\)，那有没有什么能够代替单位复根且解决精度问题呢？这就是原根。

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

若\(p\)为质数，设\(g\)为\(p\)的原根，那么\(g^i \mod p(1<j<p,1\le i\le 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,1\le i\le p-1)\)的结果两两不同。

性质

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

\(\omega_n^t\)互不相同，保证点值的合法性；
\(\omega_{2n}^{2k} = \omega_n^k\)，用于分治；
\(\omega_n^{k+\frac{n}{2}} = -\omega_n^k\)，用于分治；
当\(k\neq 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} \times 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}}\neq \omega_n^0\)，即\(\omega_n^{\frac{n}{2}}\equiv -1(\mod p)\)。可推出\(\omega_n^{k+\frac{n}{2}}=\omega_n^k \times \omega_n^{\frac{n}{2}}=-\omega_n^k (\mod p)\)，满足 性质三 。

性质四

当\(k\neq 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-1\equiv \omega_n^n-1\equiv 0 (\mod p)\)，故\(S(\omega_n^k)=0\)，性质四成立。

求原根

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

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

代码实现

将\(FFT\)里所有关于\(\omega_n\)的运算替换成\(g^q\)在模意义下的运算即可，注意\(\div n\)要改为\(\times 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;

}