快速数论变换 NTT

模板题：P3803 【模板】多项式乘法 (FFT)

建议配合快速傅里叶变换 FFT 和阶与原根两篇文章使用。

在前两篇我们介绍了 FFT 算法如何利用单位根 \(\omega_n^k\) 的性质以快速计算 DFT ，以及原根的性质。

提示：鉴于 FFT 算法的特点，以下默认 \(n\) 为 \(2\) 的幂，即 \(n = 2^m\,(m \in \mathbb{N})\) 。

FFT 那篇可能写的比较复杂，总结一下可以发现，无非就是以下几个重要性质。

\(\omega_n^k\,(0\le k \le n-1)\) 互不相同。
\(\omega_n^{k+\frac{n}{2}} = -\omega_n^k\) 。
\(1+\omega_n^k+\left(\omega_n^k\right)^2+\cdots+\left(\omega_n^k\right)^{n-1}=0\) 。

接着我们在原根的诸多性质中找出几个特别的

对于质数 \(p=qn+1\,(n=2^m)\) ，原根 \(g\) 满足 \(g^{qn} \equiv 1 \pmod{p}\) 。
对于质数 \(p\) ，设 \(g\) 为 \(p\) 的原根，则 \(g^k \bmod p\,(0 < k < p)\) 的结果两两不同。
由第二条结论可知，\(g^{k/2}\equiv -1 \pmod{p}\) 。

这里简单证明一下第二条结论。

证明考虑反证法。设 \(g^{i}\equiv g^j \pmod{p}\,(0 < i < j < p)\) ，则 \(g^{j-i}(g^i-1) \equiv 0 \pmod{p}\)

由此可知 \(g^i \equiv 0 \pmod{p}\) ，根据定义 \(i\) 只可能等于 \(p-1\) ，而此时不存在合法的 \(j\) ，与假设矛盾。 \(\square\)

于是我们发现，原根满足作为 FFT 算法核心的条件，即 \[ \omega_n \equiv g^{(p-1) / n} \pmod{p} \] 这里 \(p\) 的定义和上面的一样。

那么我们只需要把所有的 \(\omega_n\) 全部换成 \(g^{(p-1)/n}\) 就好了，这便是快速数论变换 NTT (FNTT) 。

最后，因为任意 \(n_0 \le n\) 均可以以 \(p\) 作为模数，所以我们一般取 \(998244353\) 作为 NTT 模数。

那么，我们为什么要用 NTT 呢？因为 NTT 可以在整数范围内计算，较好的免去了精度问题。

然而普通的 NTT 仅能解决特定模数（NTT 模数）的多项式乘法，或者值域有限的多项式乘法。

时间复杂度 \(\mathcal{O}(n \log n)\)

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(1e7 + 15))

const int P = 998244353;
const int G = 3, Gi = 332748118;
int a[N], b[N];
int l, r[N], limit = 1;
int qpow(int a, int b)
{
    int r = 1;
    while(b)
    {
        if(b & 1) r = r * a % P;
        b >>= 1; a = a * a % P;
    }
    return r;
}
void NTT(int *A, int type)
{
    for(int i = 0; i < limit; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    for(int mid = 1; mid < limit; mid *= 2)
    {
        int Wn = qpow((type == 1 ? G : Gi), (P - 1) / (mid * 2));
        for(int j = 0; j < limit; j += (mid * 2))
        {
            int w = 1;
            for(int k = 0; k < mid; k++, w = (w * Wn) % P)
            {
                int x = A[j + k], y = w * A[j + k + mid] % P;
                A[j + k] = (x + y) % P; A[j + k + mid] = (x - y + P) % P;
            }
        }
    }
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n, m; cin >> n >> m;
    for(int i = 0, x; i <= n; i++) cin >> x, a[i] = (x + P) % P;
    for(int i = 0, x; i <= m; i++) cin >> x, b[i] = (x + P) % P;
    for(; limit <= n + m; l++, limit *= 2);
    for(int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    NTT(a, 1); NTT(b, 1);
    for(int i = 0; i < limit; i++) a[i] = a[i] * b[i] % P;
    NTT(a, -1);
    int Inv = qpow(limit, P - 2);
    for(int i = 0; i <= n + m; i++)
        cout << a[i] * Inv % P << " \n"[i == n + m];
    return 0;
}