从傅里叶变换(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\)，由费马小定理可知
\(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^{\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\)时
由性质三中的推论可知，\((\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;
}
