FFT(快速傅里叶变换)的介绍及其C++实现_多项式乘法傅立叶变换c++

1、FFT的介绍

快速傅里叶变换 (Fast Fourier Transform)**，即利用计算机计算离散傅里叶变换（DFT）的高效、快速计算方法的统称，简称FFT，于1965年由J.W.库利和T.W.图基提出。

对于多项式

定义其乘积fg为

显然我们可以以 O(n2) 的复杂度计算这个乘积的每一项的系数。

但FFT可以以 O(nlog⁡n) 的时间复杂度来计算这个乘积。

自1965年提出至今，FFT算法在一维的信号处理和二维的图像处理中都起到了重要的作用，在人工智能、自动化等领域都有着不可替代的作用

2、FFT的C++实现

FFT的具体实现方法有很多，这里举一个递归实现的方法

#include <bits/stdc++.h>
const int MAX_COEFFICIENTS = 1 << 22;
const double EPSILON = 1e-6, PI = acos(-1.0);
​
using namespace std;
​
complex<double> coefficientsA[MAX_COEFFICIENTS], coefficientsB[MAX_COEFFICIENTS];
​
int degreeA, degreeB;
​
// FFT: Fast Fourier Transform
// Parameters:
//   - coefficients: Array representing the coefficients of a polynomial
//   - degree: Degree of the polynomial
//   - inv: If inv is 1, perform FFT; if inv is -1, perform IFFT
void FFT(complex<double> *coefficients, int degree, int inv) {
    if (degree == 1) // If there is only one term, return directly
        return;
​
    int mid = degree / 2;
​
    // Split the polynomial into two parts
    complex<double> A1[mid + 1], A2[mid + 1];
    for (int i = 0; i <= degree; i += 2) {
        A1[i / 2] = coefficients[i];
        A2[i / 2] = coefficients[i + 1];
    }
​
    // Recursively apply FFT or IFFT
    FFT(A1, mid, inv);
    FFT(A2, mid, inv);
​
    complex<double> w0(1, 0), wn(cos(2 * PI / degree), inv * sin(2 * PI / degree)); // Unit root
​
    // Merge the polynomials
    for (int i = 0; i < mid; ++i, w0 *= wn) {
        coefficients[i] = A1[i] + w0 * A2[i];
        coefficients[i + degree / 2] = A1[i] - w0 * A2[i];
    }
}
​
int main() {
    scanf("%d%d", &degreeA, &degreeB);
​
    // Input coefficients for the first polynomial
    for (int i = 0; i <= degreeA; ++i) {
        double x;
        scanf("%lf", &x);
        coefficientsA[i].real(x);
    }
​
    // Input coefficients for the second polynomial
    for (int i = 0; i <= degreeB; ++i) {
        double x;
        scanf("%lf", &x);
        coefficientsB[i].real(x);
    }
​
    int length = 1 << max((int)ceil(log2(degreeA + degreeB)), 1); // Determine the length for FFT
    FFT(coefficientsA, length, 1); // Convert coefficient representation to point-value representation
    FFT(coefficientsB, length, 1);
​
    // Perform polynomial multiplication in O(n) time
    for (int i = 0; i <= length; ++i)
        coefficientsA[i] = coefficientsA[i] * coefficientsB[i];
​
    FFT(coefficientsA, length, -1); // Convert point-value representation back to coefficient representation
​
    // Output the result
    for (int i = 0; i <= degreeA + degreeB; ++i)
        printf("%.0f ", coefficientsA[i].real() / length + EPSILON); // Divide by length; EPSILON to handle precision issues
​
    return 0;
}
​

以下是针对这个递归实现的优化

#include <bits/stdc++.h>
const int MAX_COEFFICIENTS = 1 << 22;
const double EPSILON = 1e-6, PI = acos(-1.0);
​
using namespace std;
​
complex<double> coefficientsA[MAX_COEFFICIENTS], coefficientsB[MAX_COEFFICIENTS];
​
int degreeA, degreeB;
​
// FFT: 快速傅里叶变换
// 参数：
//   - coefficients: 表示多项式的系数数组
//   - degree: 多项式的阶数
//   - inv: 如果 inv 为 1，执行 FFT；如果 inv 为 -1，执行 IFFT
void FFT(complex<double> *coefficients, int degree, int inv) {
    // 位反转置换
    for (int i = 1, j = 0; i < degree; ++i) {
        for (int k = degree >> 1; (j ^= k) < k; k >>= 1);
        if (i < j) swap(coefficients[i], coefficients[j]);
    }
​
    // Cooley-Tukey算法迭代过程
    for (int len = 2; len <= degree; len <<= 1) {
        double ang = inv * 2 * PI / len;
        complex<double> wlen(cos(ang), sin(ang));
​
        for (int i = 0; i < degree; i += len) {
            complex<double> w(1, 0);
​
            for (int j = 0, mid = len >> 1; j < mid; ++j) {
                complex<double> u = coefficients[i + j], v = w * coefficients[i + j + mid];
                coefficients[i + j] = u + v;
                coefficients[i + j + mid] = u - v;
                w *= wlen;
            }
        }
    }
}
​
int main() {
    scanf("%d%d", &degreeA, &degreeB);
​
    // 输入第一个多项式的系数
    for (int i = 0; i <= degreeA; ++i) {
        double x;
        scanf("%lf", &x);
        coefficientsA[i].real(x);
    }
​
    // 输入第二个多项式的系数
    for (int i = 0; i <= degreeB; ++i) {
        double x;
        scanf("%lf", &x);
        coefficientsB[i].real(x);
    }
​
    int length = 1 << max((int)ceil(log2(degreeA + degreeB)), 1); // 确定 FFT 的长度
    FFT(coefficientsA, length, 1); // 将系数表示转换为点值表示
    FFT(coefficientsB, length, 1);
​
    // 在 O(n) 时间内执行多项式乘法
    for (int i = 0; i <= length; ++i)
        coefficientsA[i] = coefficientsA[i] * coefficientsB[i];
​
    FFT(coefficientsA, length, -1); // 将点值表示转换回系数表示
​
    // 输出结果
    for (int i = 0; i <= degreeA + degreeB; ++i)
        printf("%.0f ", coefficientsA[i].real() / length + EPSILON); // 除以长度；EPSILON 用于处理精度问题
​
    return 0;
}
​

优化点：

1. 位反转置换（Bit-reverse permutation）： 避免了递归过程中的重复计算，通过位反转置换一次性完成了重新排列。

2. Cooley-Tukey算法的迭代实现： 使用迭代而非递归方式，减少了递归调用的开销。

3. 减少内存使用： 不再需要额外的临时数组，直接在原数组上进行计算，减少了内存占用。