多项式乘法的迭代快速傅立叶变换

给定两个多项式 A(x) 和 B(x)，求乘积 C(x) = A(x)*B(x)。在上一篇文章中，我们讨论了递归方法来解决这个复杂度为 O(nlogn) 的问题。
例子：

输入：A[] = {9, -10, 7, 6} B[] = {-5, 4, 0, -2} 输出：C(x) = $-12x^6 -14x^5 +44x^4 -20x^3 - 75x^2 +86x -45$

在信号处理等现实生活应用中，速度非常重要，本文研究了一种高效的 FFT 实现。本文重点介绍 FFT 算法的迭代版本，该算法在 O(nlogn) 时间内运行，但可以具有比递归版本更低的隐藏常数，并且可以节省递归堆栈空间。
先决条件：FFT的递归算法。
回忆上一篇文章中的递归 FFT 伪代码，在 for 循环评估中$y_k$ , $w_n^k$计算两次。我们可以将循环更改为只计算一次，将其存储在临时变量 t 中。所以，它变成，
叉$\leftarrow$ 0 到 n/2 – 1
做$\leftarrow \omega y^{\left [ 1 \right ]}$
$y_{k}\leftarrow \omega y_{k}^{\left [ 0 \right ]} + t$
$y_{k + \left ( n/2 \right )}\leftarrow _{k}^{\left [ 0 \right ]} - t$
$\omega \leftarrow \omega \omega _{n}$
此循环中的操作，乘以旋转因子 w = $w_n^k$经过$y_k^[^1^]$ ，将乘积存入 t 中，并从中加减 t $y_k^[^0^]$ , 称为蝶式运算。
从图形上看，这就是蝴蝶操作的样子：

让我们取 n=8 并继续形成迭代 fft 算法。查看上面的递归树，我们发现如果我们将初始系数向量 a 的元素按照它们出现在叶子中的顺序排列，我们可以跟踪递归 FFT 过程的执行，但是自底向上而不是顶部下。首先，我们取元素对，使用一个蝶形运算计算每一对的 DFT，然后用它的 DFT 替换这对元素。然后该向量包含 n/2 个 2 元素 DFT。接下来，我们成对获取这些 n/2 个 DFT，并通过执行两个蝶形运算来计算它们来自的四个向量元素的 DFT，用一个 4 元素 DFT 替换两个 2 元素 DFT。然后该向量包含 n/4 个 4 元素 DFT。我们以这种方式继续，直到向量包含两个 (n/2) 元素 DFT，我们使用 n/2 蝶形运算将它们组合成最终的 n 元素 DFT。

为了将这种自下而上的方法转化为代码，我们使用了一个数组 A[0…n]，它最初按照它们在树的叶子中出现的顺序保存输入向量 a 的元素。因为我们必须在树的每个级别上组合 DFT，所以我们引入了变量 s 来计算级别，范围从 1(在底部，当我们组合成对以形成 2 元素 DFT 时)到 lgn(在顶部，当佩戴时结合两个 n/2 元素 DFT 以产生最终结果)。因此算法是：

1. s=1 到 lgn 2. k=0 到 n-1 $2^s$ 3. 将两者结合起来$2^s-1$ -A[k…k+ 中的元素 DFT $2^s-1$ -1] 和 A[k+ $2^s-1$ …k+ $2^s$ -1] 到 A[k…k+ 中的一个 2s 元素 DFT $2^s$ -1]

现在为了生成代码，我们按照叶子的顺序排列系数向量元素。示例 – 叶子 0, 4, 2, 6, 1, 5, 3, 7 的顺序是索引的一点反转。从000, 001, 010, 011, 100, 101, 110, 111开始，将每个二进制数的位取反得到000, 100, 010, 110, 001, 101, 011, 111.伪码为FFT：

BIT-REVERSE-COPY(a, A) n = length [a] for k = 0 to n-1 do A[rev(k)] = a[k] ITERATIVE-FFT BIT-REVERSE-COPY(a, A) n = length(a) for s = 1 to log n do m= $2^s$$w_m$ = $e^(2*PI*i/m)$对于 j = 0 到 m/2-1 为 k = j 到 n-1 由 m 做 t = $w$ A[k+m/2] u = A[k] A[k] = u+t A[k+m/2] = ut $w = w*w_n$返回 A

从下面的并行 FFT 电路中会更清楚：

CPP

// CPP program to implement iterative
// fast Fourier transform.
#include 
using namespace std;
 
typedef complex cd;
const double PI = 3.1415926536;
 
// Utility function for reversing the bits
// of given index x
unsigned int bitReverse(unsigned int x, int log2n)
{
    int n = 0;
    for (int i = 0; i < log2n; i++)
    {
        n <<= 1;
        n |= (x & 1);
        x >>= 1;
    }
    return n;
}
 
// Iterative FFT function to compute the DFT
// of given coefficient vector
void fft(vector& a, vector& A, int log2n)
{
    int n = 4;
 
    // bit reversal of the given array
    for (unsigned int i = 0; i < n; ++i) {
        int rev = bitReverse(i, log2n);
        A[i] = a[rev];
    }
 
    // j is iota
    const complex J(0, 1);
    for (int s = 1; s <= log2n; ++s) {
        int m = 1 << s; // 2 power s
        int m2 = m >> 1; // m2 = m/2 -1
        cd w(1, 0);
 
        // principle root of nth complex
        // root of unity.
        cd wm = exp(J * (PI / m2));
        for (int j = 0; j < m2; ++j) {
            for (int k = j; k < n; k += m) {
 
                // t = twiddle factor
                cd t = w * A[k + m2];
                cd u = A[k];
 
                // similar calculating y[k]
                A[k] = u + t;
 
                // similar calculating y[k+n/2]
                A[k + m2] = u - t;
            }
            w *= wm;
        }
    }
}
 
int main()
{
    vector a{ 1, 2, 3, 4 };
    vector A(4);
    fft(a, A, 2);
    for (int i = 0; i < 4; ++i)
        cout << A[i] << "\n";
}

Input:  1 2 3 4
Output:
(10, 0)
(-2, -2)
(-2, 0)
(-2, 2)

时间复杂度分析：
复杂度是 O(nlgn)。为了说明这一点，我们将 O(nlgn) 时间内运行的最内层循环显示为：
$L(n) = \sum_{s=1}^{lg \ n}\sum_{j=0}^{2^{s-1}-1}\frac{n}{2^{s}}$
$=\sum_{s=1}^{lg \ n}\frac{n}{2^{s}}.2^{s-1}$
$=\sum_{s=1}^{lg \ n}\frac{n}{2}$
$=\Theta (n \ lg \ n)$

如果您希望与专家一起参加现场课程，请参阅DSA 现场工作专业课程和学生竞争性编程现场课程。