快快TLE

pan_g · 2024-02-16 21:04:08 · 算法·理论

快速傅里叶变换 (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)};
    // 此处不加负号是因为本质上是一样的，只用保证 IFFT 的时候加负号就好了，其实 1/n 的系数也可以都转化为 1/{\sqrt n}
    // 求出的单位根
    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), FFT(R, len>>1);
    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], [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 怎么求：

其他思路由于均来自 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 了。