理论知识

设 g 是质数 p 的一个 原根，即 1, g, g^2, \dots, g^{p - 2} 在模 p 意义下互不相同，且 g^{p - 1} \equiv 1 \pmod{p}。

将 p - 1 分解为 2 的幂与其他数的乘积，有 p - 1 = 2^l \cdot r。令 \alpha_N = g^{r}，N = 2^l，则 \alpha 有着与复数单位根 \omega_N 非常相似的性质：

而应用 FFT 加速多项式乘法，可以看做向多项式代入了 n 个点，转化为点值表示法。NTT 同理，虽然 NTT 变换后得到的“频谱”不再具有物理意义，但依然可以视为在模 p 意义下向多项式代入了 n 个点。

算法

描述

对于质数 p 和原根 g，保证 p - 1 是 n = 2^l 的倍数的情况下，只需将 FFT 中的 \omega = \mathrm{e}^{\mathrm{i} \frac{2 \pi}{n}} 替换为 g^{\frac{p - 1}{n}} 即可。

常用模数：

2^{21} > 2 \times 10^6，一般数组长度不会超过 1 \times 10^6，因此这些模数已经足够。

其余实现与 快速 Fourier 变换 几乎完全相同。

代码

ull tmp[MAXN], f[MAXN], g[MAXN], h[MAXN], omega[MAXN + 1];

void init()
{
    omega[0] = 1;
    constexpr size_t n = MAXN;
    ull o = qpow(3, (P - 1) / n);
    for (size_t i = 1; i <= n; i++)
    {
        omega[i] = omega[i - 1] * o % P;
    }
}

void ntt(ull *f, size_t n, int rev)
{
    if (n == 1)
    {
        return;
    }
    for (size_t i = 0; i < n; i++)
    {
        tmp[i] = f[i];
    }
    for (size_t i = 0; i < n / 2; i++)
    {
        f[i] = tmp[i * 2]; // g
        f[i + n / 2] = tmp[i * 2 + 1]; // h
    }
    ntt(f, n / 2, rev);
    ntt(f + n / 2, n / 2, rev);
    for (size_t i = 0; i < n / 2; i++)
    {
        ull o = omega[MAXN * (rev == -1) + rev * i * (MAXN / n)];
        tmp[i] = (f[i] + o * f[i + n / 2]) % P;
        tmp[i + n / 2] = (P + f[i] - o * f[i + n / 2] % P) % P;
    }
    for (size_t i = 0; i < n; i++)
    {
        f[i] = tmp[i];
    }
}

多模数 NTT

有些题目的数据范围或者答案会超过常用模数的范围，例如达到 10^{18}，不能直接应用 NTT。此时，可以使用多个模数同时计算 NTT。最常用的就是两个模数。

双模数 NTT 下，数据对两个模数取模后，分别对两个模数使用 NTT 计算卷积。计算完成后，得到的答案是对两个模数分别取模的结果，使用 中国剩余定理 就可以从两个取模数还原出答案。

设答案对 p_1, p_2 取模的结果为 a_1, a_2，则要寻找原来的答案 x，就是解同余方程组：

由 CRT 的公式，还原得到的数为

其中 p_2^{-1} = p_2^{p_1 - 2} \bmod p_1，p_1^{-1} = p_1^{p_2 - 2} \bmod p_2。