CC++ & Algorithm

FFT快速傅里叶变换

较难9
语言版本:C++
概述:用“果汁配方”的故事教你如何快速把两种味道的果汁混合起来,再分开找回原来的味道。

果汁配方与快速乘法:用FFT把多项式“混合”再“分离”

你有没有试过把草莓汁和蓝莓汁混在一起,变成一杯超好喝的“蓝草莓汁”?几小时后,想从混合果汁里找回原来的两种味道,该怎么办?在数学里,多项式就像不同口味的果汁,而它们的乘法就像是把果汁混合。FFT(快速傅里叶变换)就是一台神奇的“果汁分离机” —— 它能快速把两个多项式乘起来,也能快速从乘积里找回原来的多项式。

在实际编程中,FFT最常用的场景是快速计算两个多项式的乘积。比如你要计算 (1 + 2x) * (3 + 4x),普通方法要算4次乘法(1×3, 1×4, 2×3, 2×4),而FFT只用 O(n log n) 步就能搞定,当多项式次数很高时,这比普通乘法快几十倍。


1. 果汁配方:多项式其实是“味道系数”

把多项式想象成果汁配方:

  • 草莓汁:1 + 2x
    意思是:1杯纯草莓 + 2勺草莓精(x代表“勺浓度”)
  • 蓝莓汁:3 + 4x
    意思是:3杯纯蓝莓 + 4勺蓝莓精

混合后的“蓝草莓汁”配方应该是多少?答案就是两个多项式相乘的结果:

(1 + 2x) * (3 + 4x) = 1×3 + 1×4x + 2x×3 + 2x×4x
                     = 3 + 4x + 6x + 8x²
                     = 3 + 10x + 8x²

所以混合果汁的配方是:3杯纯果汁,10勺混合精,8勺浓缩双倍精(x²)。
FFT就是帮我们用更少步骤算出这个配方的“搅拌机”。


2. 普通乘法 vs. FFT 乘法:谁更快?

假设你有一大堆果汁配方(多项式有几千项),比如要计算:

(1 + 2x + 3x² + ... + 1000x⁹⁹⁹) * (2 + 4x + 6x² + ... + 2000x⁹⁹⁹)
  • 普通乘法:每个项都要配对乘一次,总共约 1000×1000 = 1,000,000 次乘法。
  • FFT乘法:先把多项式变成一堆“点值”(好比每种浓度的口味得分),然后点乘(配对乘),再变回系数。只需要大约 1000×log₂(1000) ≈ 1000×10 = 10,000 次操作。

FFT的原理用一句话说就是:把多项式从“系数表示”转换到“点值表示”,在点值空间里乘法很简单(对应位置相乘),然后快速转回系数。而这个“转来转去”的过程,就用到了分治技巧。


3. 分治思想:像分果汁一样拆开问题

FFT的核心是一个分治算法。我们拿一个多项式 A(x) = a₀ + a₁x + a₂x² + a₃x³ 来举例。

  1. 把偶数次项和奇数次项分开:

    • 偶数部分:Aeven(x) = a₀ + a₂x
    • 奇数部分:Aodd(x) = a₁ + a₃x
  2. 利用一个数学技巧:
    A(x) = Aeven(x²) + x·Aodd(x²)

  3. 对于单位根(特别选的点)可以递归计算,使得每次子问题规模减半。

这就是代码里 for (int len = 2; len <= n; len <<= 1) 所做的事情:从长度为2的小段开始,逐步合并。


4. 代码解读:一步一步看懂“果汁分离机”

下面这段C++代码实现了FFT的迭代版本(非递归),并用它计算两个多项式的乘积。每一行变量定义都加了中文注释,方便你边读边理解。

#include <iostream>
#include <complex>   // 复数类型
#include <vector>    // 动态数组
#include <cmath>     // 数学函数
using namespace std;

typedef complex<double> cd;  // 复数类型别名,double精度
const double PI = acos(-1);  // 圆周率π

// 快速傅里叶变换(FFT)
// a: 输入的多项式系数(用复数表示),长度必须是2的幂
// invert: false表示正向变换(系数→点值),true表示反向变换(点值→系数)
void fft(vector<cd>& a, bool invert) {
    int n = a.size();   // 长度

    // 第一步:将系数重新排列(位逆序置换),相当于分治的“拆开”
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;   // 从最高位开始
        for (; j & bit; bit >>= 1) j ^= bit;  // 找下一个逆序位置
        j ^= bit;
        if (i < j) swap(a[i], a[j]);  // 交换位置
    }

    // 第二步:迭代合并,len是当前合并的小段长度(2,4,8,...)
    for (int len = 2; len <= n; len <<= 1) {
        double ang = 2 * PI / len * (invert ? -1 : 1);  // 旋转因子的角度
        cd wlen(cos(ang), sin(ang));  // 长度为len的旋转因子
        for (int i = 0; i < n; i += len) {  // 每段len个元素
            cd w(1);  // 当前旋转因子,初始为1
            for (int j = 0; j < len/2; j++) {
                cd u = a[i+j];                   // 左半部分
                cd v = a[i+j+len/2] * w;         // 右半部分乘旋转因子
                a[i+j] = u + v;                  // 合并后的左半
                a[i+j+len/2] = u - v;            // 合并后的右半
                w *= wlen;                       // 更新旋转因子
            }
        }
    }

    // 如果是反向变换,要把每个数除以n(恢复原来的系数)
    if (invert) for (cd &x : a) x /= n;
}

// 两个整数系数的多项式乘法,返回结果系数向量
vector<int> multiply(vector<int>& a, vector<int>& b) {
    // 把整数系数转为复数系数
    vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    int n = 1;
    // 找到大于等于a+b长度的最小2的幂(FFT要求长度是2的幂)
    while (n < a.size() + b.size()) n <<= 1;
    fa.resize(n);   // 补0到长度n
    fb.resize(n);

    // 正向变换(系数→点值)
    fft(fa, false);
    fft(fb, false);

    // 点值乘法:对应位置相乘
    for (int i = 0; i < n; i++) fa[i] *= fb[i];

    // 反向变换(点值→系数)
    fft(fa, true);

    // 将复数结果取整并转成整数
    vector<int> result(n);
    for (int i = 0; i < n; i++) result[i] = round(fa[i].real());
    return result;
}

int main() {
    vector<int> a = {1, 2}; // 表示多项式 1 + 2x
    vector<int> b = {3, 4}; // 表示多项式 3 + 4x
    auto res = multiply(a, b);  // 调用乘法
    // 输出结果:3 10 8 对应 3 + 10x + 8x²
    for (int x : res) cout << x << " ";
    return 0;
}

运行结果3 10 8 0(最后一个0是补位,可以忽略)。


5. 常见错误与避坑指南

新手在使用FFT时容易掉进下面几个坑:

坑1:长度不是2的幂

FFT要求输入长度必须是2的幂次。比如多项式次数为3(共4项),下一个2的幂是4,刚好;但如果次数为4(共5项),就要补到8。代码中用 while (n < a.size() + b.size()) n <<= 1; 自动算出了最接近的2的幂。

坑2:忘记除以n

反向变换(invert=true)后,每个点值要除以n才能还原回原始的系数。如果漏掉,结果会放大n倍。代码里 if (invert) for (cd &x : a) x /= n; 就是做这件事。

坑3:浮点误差

FFT用double计算,乘法结果会有微小误差,比如本应是3.0的结果变成2.9999999。用 round() 取整可以修正,但如果系数很大或精度不够,可能出错。更稳妥的做法是用 long double 或优化NTT(数论变换)来避免浮点。

坑4:忘记处理负数系数

如果多项式有负系数,结果也可能出现负数,取整时需要用 round 而不是 floorround() 会四舍五入,适合处理浮点接近整数的情况。


6. 完整示例:计算 (1 + 2x + 3x²) × (4 + 5x)

我们用代码算一算,并手动验证。

int main() {
    vector<int> a = {1, 2, 3}; // 1 + 2x + 3x²
    vector<int> b = {4, 5};    // 4 + 5x
    auto res = multiply(a, b);
    // 期望结果:4*1=4, 1*5+2*4=5+8=13, 2*5+3*4=10+12=22, 3*5=15
    // 即:4 + 13x + 22x² + 15x³
    for (int x : res) cout << x << " ";  // 输出:4 13 22 15 0
    return 0;
}

输出 4 13 22 15 0,和手动计算一致(末尾0是补位)。


7. 相关知识点指引

如果你想继续深入学习,可以看看这些:

  • 复数与单位根:FFT为什么用复数的单位根?因为单位根有很好的性质(比如旋转、共轭),能让分治顺利进行。
  • 分治与递归:FFT的分治思路和归并排序很像,理解归并排序后更容易掌握。
  • NTT(数论变换):用整数模的根代替复数,避免浮点误差,适合竞赛中需要精确整数的场景。
  • 卷积:FFT本质是快速计算卷积(convolution),多项式乘法就是卷积的一种。在信号处理、图像模糊、字符串匹配中都有应用。

FFT是算法竞赛中的“重型武器”,虽然实现稍复杂,但一旦掌握,就能轻松处理上万次多项式的乘法,像变魔术一样快速“混合”又“分离”果汁配方。多写几个例子,试试不同的多项式,你就能慢慢玩转它了。

例题精讲

1单选题

FFT算法的时间复杂度是?

AO(n)
BO(n log n)
CO(n^2)
DO(log n)
2判断题

Cooley-Tukey FFT算法要求输入序列长度为2的幂。

3填空题
在递归FFT实现中,补全旋转因子:complex<double> w = ___;