快速傅里叶变换 (Fast Fourier Transform)

Fast Fast TLE

By pan_g

Update on 2025/05/10 ： 感觉之前写的比较烂，稍微改了一下，变得不那么通俗易懂了（其实还好

Update on 2025/07/15 ： 加了一点小巧思，和 MTT 。

前置芝士

包括复数、多项式、矩阵等。

离散傅里叶变换

简称 DFT ，是一种常用的简化计算手段。

对于一个有特殊性质的有限数列 $ {x_n}$ 可以考虑采用 DFT 的方式简化计算（一般在 MO 中使用的多一些吧）。

形式为

后面是 $n$ 次单位根的形式，相当于就是把 $ {x_n} $ 的生成函数，把单位根 $ e^{-i \frac{2\pi}{n}k} $ 带入的点值存入 $ X_k $ 。

对于其逆变换 IDFT ，形式为

其实是很像的，对于连续的傅里叶变换，就可以改成积分形式，并将 $n$ 设为周期即可。

多项式

对于两个多项式

便有 $h(x) = f(x)g(x) = 1x^6+10x^5+14x^4+51x^3+58x^2+41x^1+45x^0$

FFT 干的事情本质上是将多项式转化为点值表达，可用于求多项式的加卷积，上述操作记为 $h = f \circ g$ ，也有记为 $h = f \oplus g$ 的。

也可以观察到把 $10$ 带进去就是两个大数相乘的值，所以 FFT 也有优化高精度乘法的作用。

引入

人显然是可以拥有惊人的注意力的，所以

可以观察到，离散傅里叶变换的形式，看起来就是一个矩阵：

记 $\omega_n = e^{-i\frac{2\pi}{n}}$

但是，这个方法并不高效，甚至是 $O(n ^ 3)$ 的算法。

FFT

我们可以思考一些莫名其妙的性质：

1. $\omega_{\frac{n}{2}}^k = \omega_n^{2k}$
2. $\omega_n^n = 1$
3. $\omega_n^k = -\omega_n^{k+\frac{n}{2}}$

现有一个多项式 $A(x) = a_0x^0+a_1x^1+\cdots +a_{n-1}x^{n-1}$

举个例子，当 $n = 8$ 时，

将 $\omega_n^k$ 代入，

利用一点点的性质以及 $L(x^2)$ 与 $R(x^2)$ 为偶函数，所以：

欸，这不就做到分治了吗？于是就可以递归实现一个代码。

void FFT(Complex *a, int len){
    // len 存的是一半的长度
    if(!len) return ;
    Complex L[N], R[N];
    for(int i = 0;i < len;i++){
        L[i] = a[i<<1];
        R[i] = a[i<<1|1];
    }
    FFT(L, len>>1), FFT(R, len>>1);
    Complex Wn = (Complex){cos(Pi/len), sin(Pi/len)};
    Complex W = (Complex){1, 0};
    for(int i = 0;i < len;i++, W = W*Wn){
        a[i] = L[i]+W*R[i];
        a[i+len] = L[i]-W*R[i];
    }
    return ;
}

注意一点，十分重要，就是 FFT 中多项式的项数必然是 $2^k$ ，因为必须保证 $L(x)$ 和 $R(x)$ 大小都是强制对应的，所以高位要强制补足。

IFFT

人总是有惊人的注意力的，所以

IDFT 也可以类似的用矩阵表达

令 $\varepsilon_n = \dfrac{1}{\omega_n} = e^{i\frac{2\pi}{n}}$ ，然后可以观察到：

可以看出这是单位矩阵的 $n$ 倍，接着可以得到，

跟上面一样，可以放心大胆地直接复用 FFT 求逆。

void FFT(Complex *a, int len, int inverse){
    // len 存的是一半的长度
    if(!len) return ;
    Complex L[N], R[N];
    for(int i = 0;i < len;i++){
        L[i] = a[i<<1];
        R[i] = a[i<<1|1];
    }
    FFT(L, len>>1, inverse), FFT(R, len>>1, inverse);
    Complex Wn = (Complex){cos(Pi/len), inverse*sin(Pi/len)};
    Complex W = (Complex){1, 0};
    for(int i = 0;i < len;i++, W = W*Wn){
        a[i] = L[i]+W*R[i];
        a[i+len] = L[i]-W*R[i];
        if(inverse==-1){
            a[i].real/=(len<<1);
            a[i+len].real/=(len<<1);
        }
    }
    return ;
}

单位根的倒数就是它的共轭复数，这个用诱导公式推一下就好了，非常简单。

优化 FFT

因为递归导致的一些莫名其妙的耗时耗空间，所以说用迭代优化 FFT 是必然的。

对于八个系数 $a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7$ ，经过几次变换操作后，位置是这样的。

第一次： $a_0 a_2 a_4, a_6, a_1 a_3 a_5 a_7$

第二次： $a_0 a_4, a_2 a_6, a_1 a_5, a_3 a_7$

第三次： $a_0 a_4 a_2 a_6 a_1 a_5 a_3 a_7$

用逗号分隔，表示每次运算系数要计算对应的单位根的几次方。

经过观察 （观察力惊人） ，它的对应位置与对应下标的二进制正好相反，即位逆序变换。

小证一下： 我们发现对于二进制最后一位为 $1$ 的，最后排在了后 $2^{n-1}$ 位，对于其中任意一个 $2^{n-1}$ 的数列，再进行观察，就是前 $2^{n-2}$ 位，都是二进制倒数第 $2$ 位为 $1$ 的数……

相当于对于 $7$ 在 $16$ 长的数列里，第一次变换会使它在后 $8$ 位，因为它是奇数，而 $8$ 以后的数第一位都是 $1$ ，根据 $R(x)$ 的定义， $7$ 在后 $8$ 位中，地位就是 $\dfrac{7-1}{2}$ ，也就是 $3$ 或者可以理解为代码中 7 >> 1 ，也就是二进制倒数第二位而 $3$ 也是奇数，所以它就会放在后 $8$ 位的后 $4$ 位，以此类推，就可以知道， $7$ 是放在 $14$ 的位置上。

思考一下，便可得证。

同样根据证明的方式或者位逆序变换的意思，可以得出其位置的递推式。

然后可以推出来每个对应位置怎么计算。

void FFT(Complex *a, int inverse){
    for(int i = 0;i < Len;i++){ //Len是指总长度
        if(i < Reverse[i]) swap(a[i], a[Reverse[i]]);
    } //if语句是防止多次交换
    for(int L = 1;L < Len;L<<=1){ //L表示遍历到现在的长度的一半
        Complex Wn = (Complex){cos(Pi/L), sin(Pi/L)*inverse};
        for(int l = 0;l < Len;l += (L << 1)){
            Complex W = (Complex){1, 0};
            for(int p = 0;p < L;p++, W = W*Wn){
                Complex x = a[l+p], y = a[L+l+p]*W;
                a[l+p] = x+y, a[l+p+L] = x-y;
                //蝶形运算的优化（可以节省一个较大的常数）
            }
        }
    } // 迭代递推FFT
    if(inverse == -1) for(int i = 0;i < Len;i++) a[i].real /= Len;
    return ;
}

FFT 就到此结束了，最后求值，就长这样。

// 将两个多项式的系数扔进 a, b 数组中
while(Len < n+m+1) Len <<= 1, x++;
FFT(a, 1), FFT(b, 1); //FFT(a), FFT(b)
for(int i = 0;i < Len;i++) a[i] = a[i]*b[i];
FFT(a, -1); //IFFT(a)
for(int i = 0;i <= Len;i++) Ans[i] = a[i].real+0.5; //此处防止出现精度误差

当然 FFT 的常数是十分大的，所以还有更快的也会减少精度误差的 NTT ，也是基于 FFT 的算法。

NTT

或者更准确的应该叫 FNTT ，全称快速数论变换，也是 OI 中最常写的。

前置芝士

原根、逆元、 FFT 等。

引入

由于 FFT 采用的是浮点运算，时间消耗会对应的长，并且会丢失精度，所以实在难以较好地被接受，在此基础上，于是就有了 NTT 。

但是，一切原因都需归咎于带入的点值的选择上。但在整个数域上，很难到第二个这么好的数，但在有限的集合内，就会很好研究。

所以，造成了 NTT 使用的最大前提，就是必须再模意义下使用，但是是可以创造模意义。

NTT

模意义下，有一个非常棒的取值方式就是取模数 $p$ 的原根 $g$ ，其中 $p$ 为质数，因为原根 $g^k \bmod p$ 的值各不相同。

令 $p = qn+1(n = 2^m)$ ，就一定有 $g^{qn} \equiv 1 \pmod p$ ，令 $g_n \equiv g^q \pmod p$ ，并将 $g_n$ 等价于 $\omega_n$ ，只因也满足重要的一些性质：

1. $g_n^n \equiv 1 \pmod p$
2. $g_{2n}^{2k} \equiv g_n^k \pmod p$
3. $g_n^{k+\frac{n}{2}} \equiv -g_n^k \pmod p$

考虑 $g_n$ 怎么求： $$ g_n = g^q = g^{\frac{p-1}{n}} $$

其他思路由于均来自 FFT ，所以其他不变即可。

void NTT(ll *a, int inverse){
    for(int i = 0;i < Len;i++){
        if(i < Reverse[i]) swap(a[i], a[Reverse[i]]);
    }
    for(int L = 1;L < Len;L<<=1){
        ll Gn = quick_pow(inverse == 1 ? g : g_, (mod-1)/(L << 1));
        // g代表模数的原根  _g代表模数的原根在模意义下的倒数
        for(int l = 0;l < Len;l += (L << 1)){
            ll G = 1;
            for(int p = 0;p < L;p++, G = G*Gn%mod){
                ll x = a[l+p]%mod, y = a[L+l+p]*G%mod;
                a[l+p] = (x+y)%mod, a[l+p+L] = (x-y+mod)%mod;
                // 注意会出现负数
            }
        }
    }
    if(inverse == -1){
        ll inv = quick_pow(Len, mod-2);
        for(int i = 0;i < Len;i++) a[i] = a[i]*inv%mod;
    }
    return ;
}

其余同 FFT 换个函数名即可。

原根性质证明： $g_n\equiv g^{\frac{p-1}{n}} \pmod p$

显然，根据费小 $g_n^n\equiv g^{p-1} \equiv 1\pmod p$

可得， $g_n^k\equiv g^{\frac{k(p-1)}{n}}\pmod p$ 。

上下同乘以 $2$ 答案不变，所以 $g_n^k\equiv g_{2n}^{2k}\pmod p$ 。

对于第三个性质， $g_n^{k+\frac{n}{2}}\equiv -g_n^k\pmod p$

要证即证， $g^{\frac{p-1}{2}}\equiv p-1\pmod p$

同时平方， $g^{p-1}\equiv p^2-2p+1\pmod p$

根据费小，显然成立。

对于本质的小巧思

众所周知，白舞鞋的花语是忘记一切，黑舞鞋的花语是我有一计。

由于先前提到过， FFT 是将多项式转化为点值表达的过程。

所以，对于这个式子：

我们并不需要乘两次，然后再减。

我们只需要对于这个式子中的 $F(x)$ 与 $G(x)$ 先 FFT 一次，得到 $F'(x)$ 与 $G'(x)$ 。

然后 $H'(x) = F'(x)G'^2(x) - 2G'(x)$ ，对于 $H'(x)$ 再 IFFT 回去就好了。

因为点值表达时，只需要把点值计算出来即可，所以这样 IFFT 是可行的。

对于任意模数的 NTT

由于，我们可能遇到的模数可能很大，所以如果直接 FFT 的话最大的数值为 $P^2N$ ，其中 $N$ 是多项式项数， $P$ 是模数。

按 $N = 10^5, P = 10^9$ 计，则这个数最大是 $10^{23}$ ，用 long long 存不下，所以可能要用 __int128_t 这种变量类型。

因此有一种更加暴力的做法就是三模 NTT ，先用三个已知原根的质数（使用三个是保证三个模数相乘能够大于直接相乘的值域大小），做 NTT 再通过 CRT 把原先的数求出来即可。

但是这样就需要 9 次 NTT ，常数巨大，令人望而却步。

于是，我们需要干一件大事，就是把原多项式的系数拆开，即令 $F(x) = kA(x) + B(x), G(x) = kC(x) + D(x)$ ，其主要目的在于让系数的范围降下来，所以一般来讲 $k = \lfloor \sqrt P \rfloor$ ，也会取相近的 $2$ 的幂次。

那么 $F(x)G(x) = k^2A(x)C(x) + k(B(x)C(x) + A(x)D(x)) + B(x)D(x)$ 。

理论上来讲， $A,B,C,D$ 各自 FFT ，再捆成三捆再 IFFT ，可以只用 7 次，但是还是太多了。

对于 IFFT ，我们知道最终结果一定是整系数多项式（因为原来进行 FFT 的多项式都是整系数），所以虚部为 $0$ ，如果将另一个点值乘上 $i$ 之后，一起 IFFT ，那么它的结果一定是在虚部，因此 IFFT 的次数可以除以 $2$ ，变为 2 次 IFFT 。

其次是对于 FFT 次数的优化，

考虑构造两个多项式 $P(x) = F(x) - iG(x)$ 和 $Q(x) = F(x) + iG(x)$ 。

对于两个多项式分别 FFT ：

接着考虑到 $\dfrac{2(n - k)j\pi}{n} = 2j\pi - \dfrac{2kj\pi}{n}$ 。

所以 $Q'(n - k) = \sum_{j = 0}^{n - 1} (A_j\cos(\dfrac{2kj\pi}{n}) + B_j\sin(\dfrac{2kj\pi}{n})) - i(A_j\sin(\dfrac{2kj\pi}{n}) - B_j\cos(\dfrac{2kj\pi}{n}))$ 。

于是我们会非常惊喜的发现 $P'(k)$ 与 $Q'(n - k)$ 共轭，但是 $P'(0)$ 与 $Q'(0)$ 共轭。

于是，我们就可以通过 2 次 FFT 加上出色的解方程芝士求出 4 个多项式的点值表达。

于是，我们就可以 2 次 FFT + 2 次 IFFT 就可以做到 4 次 FFT + 3 次 IFFT 的效果。

然后就可以实现任意模的 NTT 了。