本文已参与「新人创作礼」活动，一起开启掘金创作之路

前言

FFT是一个快速(O(nlogn))求两个多项式乘积的算法。 在阅读本文之前，请先阅读1.本文是对1的解说和弥补，首要贡献有两点：

  • 把原文没说清的理论部分补上；
  • 经过递归，把原文的思路重新理了一遍，使理论和代码都更加好了解。

别的，本文的附录里边只放了源代码，且在文中有提到。

理论：IFFT的理论根底

只要搞懂以下5个topic，就能够了解快速傅里叶变换的代码逻辑：

1.多项式的系数表明法和点值表明法 2.n次单位根 3.DFT 4.FFT 5.IFFT

理论部分的前四点构成了FFT的理论根底。 对于这四点，1都有，且很好了解。请读者阅读原文熟悉前四点理论。

下面首要讲一下第5点（IFFT的理论根底）。

FFT本质上处理了一个乘法问题

以下一切记号和1对应，并记wi=ni=exp(i∗2i/n)w_i = \omega_n^i = exp(\mathbf{i}*2\pi i/n)，对应n次单位根。（i\mathbf{i}是虚数单位。）

记n次单位根的幂矩阵W，系数向量C，值向量V如下：

W=[Wi,j=(wi)j]n∗n=[(w0)0⋯(w0)n−1(w1)0⋯(w1)n−1⋮⋱⋮(wn−1)0⋯(wn−1)n−1]W=[W_{i,j}=(w_i)^j]_{n*n}= \left[ \begin{aligned} (w_0)^0 & \cdots & (w_0)^{n-1} \\ (w_1)^0 & \cdots & (w_1)^{n-1} \\ \vdots & \ddots & \vdots \\ (w_{n-1})^0 & \cdots & (w_{n-1})^{n-1} \\ \end{aligned} \right]
C=[a0a1⋮an−1]C = \left[ \begin{aligned} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{aligned} \right]
V=[A(x0)A(x1)⋮A(xn−1)]V = \left[ \begin{aligned} A(x_0) \\ A(x_1) \\ \vdots \\ A(x_{n-1}) \end{aligned} \right]

其间WC=VWC=V.

FFT本质上处理的问题是已知W和C，核算V的乘法问题，对这个问题给出了一个O(nlogn)O(nlogn)的办法。

IFFT能够被表明成一个类似的乘法问题

IFFT能够表明成一个类似的乘法问题，并经过相似的办法处理。这一部分参阅了2.

咱们首先求W的逆矩阵。记vi=conj(wi)=exp(−i∗2i/n)v_i = conj(w_i) = exp(-\mathbf{i}*2\pi i/n)，即wiw_i的共轭。（i\mathbf{i}是虚数单位。）

记矩阵W’为

W′=[Wi,j′=(vi)j]n∗n=[(v0)0⋯(v0)n−1(v1)0⋯(v1)n−1⋮⋱⋮(vn−1)0⋯(vn−1)n−1]W’=[W’_{i,j}=(v_i)^j]_{n*n}= \left[ \begin{aligned} (v_0)^0 & \cdots & (v_0)^{n-1} \\ (v_1)^0 & \cdots & (v_1)^{n-1} \\ \vdots & \ddots & \vdots \\ (v_{n-1})^0 & \cdots & (v_{n-1})^{n-1} \\ \end{aligned} \right]

发现WW′=W′W=nIWW’=W’W=nI（易证）,所以W的逆矩阵是W′n\frac{W’}{n},有W′V=nCW’V=nC.因而，只要在FFT的完结算法中把n次单位根改成其共轭，就能够从向量V用O(nlogn)O(nlogn)时间求解到向量C，即完结IFFT算法。

综上所述，完结IFFT和完结FFT使用的是同一套逻辑，只是使用了不同的n次单位根。

FFT算法的根底完结

以下是FFT算法的根底完结，fft_logic是FFT和IFFT的一起主逻辑，fft和ifft是给fft_logic套的两个壳。结合注释阅读fft_logic的代码即可。 （完好的程序见附录1，包含测验用的主函数。complex类的使用见3.）

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支撑位数
complex<double> b[N]; // 用来充任暂时调整空间的数组
void fft_logic(complex<double> *a, int n, int inv){
    // 参数：
    // 当inv = 1时，a是系数多项式，n是当时数组长度（2的幂次），函数作用是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当时数组长度（2的幂次），函数作用是原地变成系数多项式，可是所得的系数是n倍，需求在包装的函数中进行调整
    if (n == 1) return; // 为什么？由于omega_1^0=1，点值多项式和系数多项式的表明完全一致。
    // 使用B暂存和写回，把a的次序调整为 a[0] a[2] .. a[n-2] a[1] a[3] .. a[n-1]，前后两半
    for(int i = 0; i < n/2; i ++){
        b[i]       = a[i * 2];
        b[i + n/2] = a[i * 2 + 1];
    }
    for(int i = 0; i < n; i ++)
        a[i] = b[i];
    // 分治求A1和A2
    fft_logic(a, n/2, inv);
    fft_logic(a + n/2, n/2, inv);
    // 经过A1和A2，核算A
    double unit_rad = 2 * pi / n; // 单位角起伏值
    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i 
        complex<double> tmp1 = a[i];
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }
}
void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}
void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}

如同原文作者所说，这种办法的缺点便是常数太大需求优化…

FFT的迭代版别：优化常数

依据递归的版别

针对上一节中算法的不足，留意到ffp_logic函数中“每次递归都重排”，咱们将其优化成“预先重排后直接按次序处理”。

优化的思路是，针对原来算法中每一次A1和A2的交叉，其实能够经过提前重排原始数组来替代。拿来原文的一张图来说明：

【算法笔记】快速傅立叶变换（FFT）

n=8时，原来的次序终究被交叉成了0 4 2 6 1 5 3 7,n=4时则是0 2 1 3.这个数列仅由n决定，所以咱们能够用一个函数生成这样的交叉数列。

void fft_rearrange_decidesequence_logic(int *rev, int n){
    // 给定数组rev和数组长度n，函数的功用是将所需的次序写入数组，比方n = 4时将次序 0 2 1 3 写入数组rev
    if(n == 1){
        rev[0] = 0;
        return;
    }
    // 取得 n/2 时的次序,暂时放在rev的后半
    fft_rearrange_decidesequence_logic(rev + n/2, n/2);
    // 使用 n/2 时的次序结构 n 时的次序
    for(int i = 0; i < n/2; i ++){
        rev[i] = 2 * rev[i + n/2];
        rev[i + n/2] = 2 * rev[i + n/2] + 1;
    }
}

fft_rearrange_decidesequence_logic函数用递归的办法，自然地把交叉的逻辑描绘出来，成果得到正确的重排次序。

void fft_rearrange_logic(complex<double> *a, int n){
    // 依照算法，把a重新排列为能够直接递归向上的次序
    // 核算bit: 满意pow(2, bit) = n
    int bit = 0;
    while((1 << bit) < n) bit ++;
    // rev: 确认a终究rearrange的方位序列
    int* rev = new int[n];
    fft_rearrange_decidesequence_logic(rev, n);
    // 依照rev序列调整a的次序
    complex<double>* tmp_a = new complex<double>[n];
    for(int i = 0; i < n; i ++) tmp_a[i] = a[rev[i]];
    for(int i = 0; i < n; i ++) a[i] = tmp_a[i];
    delete[] rev;
}

fft_rearrange_logic依据n值取得重排序列，并将输入进行重排。

void fft_merge_logic(complex<double> *a, int n, int inv){
    if(n == 1) return;
    fft_merge_logic(a,       n/2, inv);
    fft_merge_logic(a + n/2, n/2, inv);
    double unit_rad = 2 * pi / n;
    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i
        complex<double> tmp1 = a[i]; 
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }
}
void fft_logic(complex<double> *a, int n, int inv){
    // 参数：
    // 当inv = 1时，a是系数多项式，n是当时数组长度（2的幂次），函数作用是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当时数组长度（2的幂次），函数作用是原地变成系数多项式，可是所得的系数是n倍，需求在包装的函数中进行调整
    // a中元素的次序调整为算法要求的次序
    fft_rearrange_logic(a, n);
    // 调整次序之后的a进行兼并
    fft_merge_logic(a, n, inv);
}

在将重排环节统一之后，fft_logic函数如上所示。在这里，fft_merge_logic把上一节中重排后兼并的逻辑抽了出来。

完好的，可运行的代码见附录2.

这一节的意图是对于原文改良速度的FFT算法写出一个易于了解，可读性强的完结。

rearrange逻辑：递归变非递归

比较一下附录2和3的代码，发现fft_rearrange_decidesequence_logic被砍掉了，换成了一个简练的循环。 （别的，依照rev序列调整a的次序的部分，尊重原文1修改了一下。我觉得两种办法效率差异不大，并且也都比较简单了解。）

由于有一个定论：每个方位分治后的终究方位为其二进制翻转后得到的方位（即老生常谈的蝴蝶算法），所以能够用这个公式直接生成rev，而不需求写递归。

void fft_rearrange_logic(complex<double> *a, int n){
    // 依照算法，把a重新排列为能够直接递归向上的次序
    // 核算bit: 满意pow(2, bit) = n
    int bit = 0;
    while((1 << bit) < n) bit ++;
    // rev: 确认a终究rearrange的方位序列
    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
    // 依照rev序列调整a的次序
    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交流两次（便是没交流）
    }
    delete[] rev;
}

完好的代码如附录3所示。

去除算法中的一切递归

咱们考虑递归的结构（宽恕我懒得再画一张新图）：

【算法笔记】快速傅立叶变换（FFT）

8依赖2个4的完结，每个4依赖其下面2个2的完结，每个2依赖下面2个1的完结。因而，能够经过一个循环，先完结8个1（不需求操作），再完结4个2，再完结2个4，最终完结8.

把这个逻辑用循环写出来，就得到了刨除递归的算法。这个算法和原文1的终究板子是对应的。

void fft_logic(complex<double> *a,int n,int inv){
    // bit: pow(2, bit) = n
    int bit = 0;
    while((1 << bit) < n) bit ++;
    // rev: 确认终究rearrange的方位序列
    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
    // 依照rev序列调整a的次序
    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交流两次（便是没交流）
    }
    delete[] rev;
    // 调整次序之后的a进行兼并
    for (int mid=1;mid<n;mid*=2){
        // mid循环内，预备把两个长度为mid的序列兼并成2 * mid的序列
        double unit_rad = 2 * pi / (2 * mid); // 单位角起伏值
        for (int i=0;i<n;i+=mid*2){
            // i循环内，把方位在i~i+mid*2方位的两个长度为mid的序列兼并成2 * mid的序列
            for (int j = 0; j < mid; j ++){
                // j循环内的逻辑和
                complex<double> x(cos(j * unit_rad), inv*sin(j * unit_rad)); // x = omega_n^i
                complex<double> tmp1 = a[i+j]; 
                complex<double> tmp2 = x * a[i+j+mid];
                a[i+j]     = tmp1 + tmp2;
                a[i+j+mid] = tmp1 - tmp2;
            }
        }
    }
}

完好的代码如附录4所示。

附录

附录1

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支撑位数
complex<double> b[N]; // 用来充任暂时调整空间的数组
void fft_logic(complex<double> *a, int n, int inv){
    // 参数：
    // 当inv = 1时，a是系数多项式，n是当时数组长度（2的幂次），函数作用是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当时数组长度（2的幂次），函数作用是原地变成系数多项式，可是所得的系数是n倍，需求在包装的函数中进行调整
    if (n == 1) return; // 为什么？由于omega_1^0=1，点值多项式和系数多项式的表明完全一致。
    // 使用B暂存和写回，把a的次序调整为 a[0] a[2] .. a[n-2] a[1] a[3] .. a[n-1]，前后两半
    for(int i = 0; i < n/2; i ++){
        b[i]       = a[i * 2];
        b[i + n/2] = a[i * 2 + 1];
    }
    for(int i = 0; i < n; i ++)
        a[i] = b[i];
    // 分治求A1和A2
    fft_logic(a, n/2, inv);
    fft_logic(a + n/2, n/2, inv);
    // 经过A1和A2，核算A
    double unit_rad = 2 * pi / n; // 单位角起伏值
    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i 
        complex<double> tmp1 = a[i];
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }
}
void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}
void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}
// 主函数测验
complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘成果
int main(){
    // 两个相乘多项式的初始化。留意：两个相乘多项式和最终的乘积多项式需求有相同的项数
    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }
    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4
    // F，乘，F逆
    fft(A1, N);
    fft(A2, N);
    /*
    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;
    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;
    */
    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];
    ifft(C, N);
    // 输出
    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;
    return 0;
}

附录2

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支撑位数
complex<double> b[N]; // 用来充任暂时调整空间的数组
void fft_merge_logic(complex<double> *a, int n, int inv){
    if(n == 1) return;
    fft_merge_logic(a,       n/2, inv);
    fft_merge_logic(a + n/2, n/2, inv);
    double unit_rad = 2 * pi / n;
    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i
        complex<double> tmp1 = a[i]; 
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }
}
void fft_rearrange_decidesequence_logic(int *rev, int n){
    // 给定数组rev和数组长度n，函数的功用是将所需的次序写入数组，比方n = 4时将次序 0 2 1 3 写入数组rev
    if(n == 1){
        rev[0] = 0;
        return;
    }
    // 取得 n/2 时的次序,暂时放在rev的后半
    fft_rearrange_decidesequence_logic(rev + n/2, n/2);
    // 使用 n/2 时的次序结构 n 时的次序
    for(int i = 0; i < n/2; i ++){
        rev[i] = 2 * rev[i + n/2];
        rev[i + n/2] = 2 * rev[i + n/2] + 1;
    }
}
void fft_rearrange_logic(complex<double> *a, int n){
    // 依照算法，把a重新排列为能够直接递归向上的次序
    // 核算bit: 满意pow(2, bit) = n
    int bit = 0;
    while((1 << bit) < n) bit ++;
    // rev: 确认a终究rearrange的方位序列
    int* rev = new int[n];
    fft_rearrange_decidesequence_logic(rev, n);
    // 依照rev序列调整a的次序
    complex<double>* tmp_a = new complex<double>[n];
    for(int i = 0; i < n; i ++) tmp_a[i] = a[rev[i]];
    for(int i = 0; i < n; i ++) a[i] = tmp_a[i];
    delete[] rev;
}
void fft_logic(complex<double> *a, int n, int inv){
    // 参数：
    // 当inv = 1时，a是系数多项式，n是当时数组长度（2的幂次），函数作用是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当时数组长度（2的幂次），函数作用是原地变成系数多项式，可是所得的系数是n倍，需求在包装的函数中进行调整
    // a中元素的次序调整为算法要求的次序
    fft_rearrange_logic(a, n);
    // 调整次序之后的a进行兼并
    fft_merge_logic(a, n, inv);
}
void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}
void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}
// 主函数测验
complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘成果
int main(){
    // 两个相乘多项式的初始化。留意：两个相乘多项式和最终的乘积多项式需求有相同的项数
    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }
    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4
    // F，乘，F逆
    fft(A1, N);
    fft(A2, N);
    /*
    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;
    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;
    */
    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];
    ifft(C, N);
    // 输出
    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;
    return 0;
}

附录3

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
const double pi = 3.14159265358979323846;
const int N = 1024; // 多项式的最大支撑位数
complex<double> b[N]; // 用来充任暂时调整空间的数组
void fft_merge_logic(complex<double> *a, int n, int inv){
    if(n == 1) return;
    fft_merge_logic(a,       n/2, inv);
    fft_merge_logic(a + n/2, n/2, inv);
    double unit_rad = 2 * pi / n;
    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i
        complex<double> tmp1 = a[i]; 
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }
}
void fft_rearrange_logic(complex<double> *a, int n){
    // 依照算法，把a重新排列为能够直接递归向上的次序
    // 核算bit: 满意pow(2, bit) = n
    int bit = 0;
    while((1 << bit) < n) bit ++;
    // rev: 确认a终究rearrange的方位序列
    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
    // 依照rev序列调整a的次序
    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交流两次（便是没交流）
    }
    delete[] rev;
}
void fft_logic(complex<double> *a, int n, int inv){
    // 参数：
    // 当inv = 1时，a是系数多项式，n是当时数组长度（2的幂次），函数作用是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当时数组长度（2的幂次），函数作用是原地变成系数多项式，可是所得的系数是n倍，需求在包装的函数中进行调整
    // a中元素的次序调整为算法要求的次序
    fft_rearrange_logic(a, n);
    // 调整次序之后的a进行兼并
    fft_merge_logic(a, n, inv);
}
void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}
void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}
// 主函数测验
complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘成果
int main(){
    // 两个相乘多项式的初始化。留意：两个相乘多项式和最终的乘积多项式需求有相同的项数
    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }
    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4
    // F，乘，F逆
    fft(A1, N);
    fft(A2, N);
    /*
    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;
    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;
    */
    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];
    ifft(C, N);
    // 输出
    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;
    return 0;
}

附录4

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支撑位数
complex<double> b[N]; // 用来充任暂时调整空间的数组
void fft_logic(complex<double> *a,int n,int inv){
    // bit: pow(2, bit) = n
    int bit = 0;
    while((1 << bit) < n) bit ++;
    // rev: 确认终究rearrange的方位序列
    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
    // 依照rev序列调整a的次序
    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交流两次（便是没交流）
    }
    delete[] rev;
    // 调整次序之后的a进行兼并
    for (int mid=1;mid<n;mid*=2){
        // mid循环内，预备把两个长度为mid的序列兼并成2 * mid的序列
        double unit_rad = 2 * pi / (2 * mid); // 单位角起伏值
        for (int i=0;i<n;i+=mid*2){
            // i循环内，把方位在i~i+mid*2方位的两个长度为mid的序列兼并成2 * mid的序列
            for (int j = 0; j < mid; j ++){
                // j循环内的逻辑和
                complex<double> x(cos(j * unit_rad), inv*sin(j * unit_rad)); // x = omega_n^i
                complex<double> tmp1 = a[i+j]; 
                complex<double> tmp2 = x * a[i+j+mid];
                a[i+j]     = tmp1 + tmp2;
                a[i+j+mid] = tmp1 - tmp2;
            }
        }
    }
}
void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}
void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}
// 主函数测验
complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘成果
int main(){
    // 两个相乘多项式的初始化。留意：两个相乘多项式和最终的乘积多项式需求有相同的项数
    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }
    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4
    // F，乘，F逆
    fft(A1, N);
    fft(A2, N);
    /*
    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;
    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;
    */
    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];
    ifft(C, N);
    // 输出
    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;
    return 0;
}

参阅

Footnotes

  1. 十分简明易懂的FFT（快速傅里叶变换） ↩ ↩2 ↩3 ↩4 ↩5 ↩6

  2. 快速傅里叶变换（FFT）和逆快速傅里叶变换（IFFT） ↩

  3. C++ complex类 ↩