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

给定两个多项式 A(x) 和 B(x)，求乘积 C(x) = A(x)*B(x)。已经有一个 O( $n^2$ ) 天真的方法来解决这个问题。这里。这种方法使用多项式的系数形式来计算乘积。
多项式的系数表示$A(x)=\sum_{j=0}^{n-1}a_jx^j$是 a = a0, a1, …, an-1。
• 例子-
$A(x) = 6x^3 + 7x^2 - 10x + 9$
$B(x) = -2x^3 + 4x - 5$
A(x) = (9, -10, 7, 6) 的系数表示
B(x) = (-5, 4, 0, -2) 的系数表示

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

如果我们以另一种形式表示多项式，我们可以做得更好。

是的

想法是以点值形式表示多项式，然后计算乘积。度数为 n 的多项式 A(x) 的点值表示是一组 n 个点值对是{ (x0, y0), (x1, y1), …, (xn-1, yn-1 )} 使得所有的 xi 都是不同的，并且 yi = A(xi) 对于 i = 0, 1, …, n-1。
例子

$A(x) = x^3 - 2x + 1$ xi — 0, 1, 2, 3 A(xi) — 1, 0, 5, 22

上述多项式的点值表示为 { (0, 1), (1, 0), (2, 5), (3, 22) }。使用霍纳的方法，(这里讨论)，n 点评估需要时间 O( $n^2$ )。它只是计算 n 个不同点的某个 x 处的 A(x) 值，因此时间复杂度为 O( $n^2$ )。现在多项式转换为点值，可以很容易地再次使用 horner 方法计算 C(x) = A(x)*B(x)。这需要 O(n) 时间。这里重要的一点是 C(x) 的度界为 2n，那么 n 个点将只给出 C(x) 的 n 个点，因此对于这种情况，我们需要 2n 个不同的 x 值来计算 2n 个不同的 y 值。现在计算出乘积，可以将答案转换回系数向量形式。为了回到系数向量形式，我们使用这个评估的逆。求值的逆过程称为插值。使用拉格朗日公式的插值将点值形式转化为多项式的系数向量形式。拉格朗日公式是 –
$A(x) = \sum^{n-1}_{k=0} y_{k} \frac{\prod _{j\neq k}(x-x_{j})}{\prod _{j\neq k}(x_{k}-x_{j})}$
到目前为止，我们讨论过，
.
这个想法还是解决了O( $n^2$ ) 时间复杂度。我们可以使用任何我们想要的点作为评估点，但是通过仔细选择评估点，我们可以在 O(n log n) 时间内在表示之间进行转换。如果我们选择“复单位根”作为评估点，我们可以通过对系数向量进行离散傅立叶变换(DFT)来产生点值表示。我们可以执行逆运算，即插值，通过对点值对进行“逆 DFT”，产生一个系数向量。快速傅立叶变换 (FFT) 可以在 O(nlogn) 时间内执行 DFT 和逆 DFT。
离散傅立叶变换
DFT 正在评估 n 个复数 n 次单位根处的多项式的值$\omega ^{0}_{n},\omega ^{1}_{n},\omega ^{2}_{n}......\omega ^{n-1}_{n}$ .因此对于$y_{k}=\omega ^{k}_{n}$ k = 0, 1, 2, …, n-1, y = (y0, y1, y2, …, yn-1) 是给定多项式的离散傅立叶变换 (DFT)。
度数为 n 的两个多项式的乘积是度数为 2n 的多项式。因此，在计算输入多项式 A 和 B 之前，我们首先通过添加 n 个为 0 的高阶系数将它们的度界加倍到 2n。因为向量有 2n 个元素，我们使用“复数 2nth 单位根”，表示为由 W2n (omega 2n)。我们假设 n 是 2 的幂；我们总是可以通过添加高阶零系数来满足这个要求。
快速傅立叶变换
这是解决这个问题的分而治之的策略。

分别使用 A(x) 的偶指数和奇指数系数，定义两个新的度数为 n/2 的多项式
$A0(x) = a0 + a2x + a4x^2 + ... + an-2x^n/2-1.$
$A1(x) = a1 + a3x + a5x^2 + ... + an-1x^n/2-1.$
$A(x) = A0(x^2) + xA1(x^2)$

评估 A(x) 的问题$\omega ^{0}_{n},\omega ^{1}_{n},\omega ^{2}_{n}......\omega ^{n-1}_{n}$简化为评估点处的度数限制 n/2 多项式 A0(x) 和 A1(x)
$(\omega ^{0}_{n})^{2},(\omega ^{1}_{n})^{2},......(\omega ^{n-1}_{n})^{2}$
现在结合结果$A(x) = A0(x^2) + xA1(x^2)$

列表$(\omega ^{0}_{n})^{2},(\omega ^{1}_{n})^{2},......(\omega ^{n-1}_{n})^{2}$不包含 n 个不同的值，但包含 n/2 个复数 n/2 个单位根。多项式 A0 和 A1 在第 n 个复数单位根处递归计算。子问题与原始问题具有完全相同的形式，但大小只有其一半。所以形成的递归是T(n) = 2T(n/2) + O(n)，即复杂度O(nlogn)。

Algorithm
1. Add n higher-order zero coefficients to A(x) and B(x)
2. Evaluate A(x) and B(x) using FFT for 2n points
3. Pointwise multiplication of point-value forms
4. Interpolate C(x) using FFT to compute inverse DFT

递归FFT的伪代码

Recursive_FFT(a){ n = length(a) // a 是输入系数向量，如果 n = 1 则返回 a // wn 是原理复数的 n 次单位根。 wn = e^(2*pi*i/n) w = 1 // 偶数索引系数 A0 = (a0, a2, …, an-2 ) // 奇数索引系数 A1 = (a1, a3, .. ., an-1 ) y0 = Recursive_FFT(A0) // 局部数组 y1 = Recursive-FFT(A1) // k = 0 到 n/2 – 1 的局部数组 // y 数组存储 DFT // of给定多项式。 do y[k] = y0[k] + w*y1[k] y[k+(n/2)] = y0[k] – w*y1[k] w = w*wn return y } 上面的递归树执行-

为什么这样做?
$y_{k} = y_{k}^{[0]} + \omega ^{k}_{n}y^{[1]}_{k}\newline y_{k} = A^{[0]}(\omega ^{2k}_{n})) + \omega ^{k}_{n}A^{[1]}(\omega ^{2k}_{n}) \newline y_{k} = A( \omega ^{k}_{n}) \newline \newline y_{k+(n/2)} = y_{k}^{[0]} - \omega ^{k}_{n}y^{[1]}_{k}\newline y_{k+(n/2)} = y_{k}^{[0]} + \omega ^{k+(n/2)}_{n} y_{k}^{[1]}\newline y_{k+(n/2)} = A^{[0]}(\omega ^{2k}_{n})) + \omega ^{k+(n/2)}_{n}A^{[1]}(\omega ^{2k}_{n})\newline y_{k+(n/2)} = A^{[0]}(\omega ^{2k+n}_{n})) + \omega ^{k+(n/2)}_{n}A^{[1]}(\omega ^{2k+n}_{n})\newline y_{k+(n/2)} = A^{[0]}(\omega ^{k+(n/2)}_{n}))$
自从， $\omega ^{k}_{n/2} = \omega ^{2k}_{n} , \omega ^{2k+n}_{n} = \omega ^{2k}_{n} , \omega ^{k+(n/2)}_{n} = -\omega ^{k}_{n}$
因此，Recursive-FFT 返回的向量 y 确实是输入的 DFT
向量

#include 
using namespace std;
  
// For storing complex values of nth roots
// of unity we use complex
typedef complex cd;
  
// Recursive function of FFT
vector fft(vector& a)
{
    int n = a.size();
  
    // if input contains just one element
    if (n == 1)
        return vector(1, a[0]);
  
    // For storing n complex nth roots of unity
    vector w(n);
    for (int i = 0; i < n; i++) {
        double alpha = -2 * M_PI * i / n;
        w[i] = cd(cos(alpha), sin(alpha));
    }
  
    vector A0(n / 2), A1(n / 2);
    for (int i = 0; i < n / 2; i++) {
  
        // even indexed coefficients
        A0[i] = a[i * 2];
  
        // odd indexed coefficients
        A1[i] = a[i * 2 + 1];
    }
  
    // Recursive call for even indexed coefficients
    vector y0 = fft(A0); 
  
    // Recursive call for odd indexed coefficients
    vector y1 = fft(A1);
  
    // for storing values of y0, y1, y2, ..., yn-1.
    vector y(n);
  
    for (int k = 0; k < n / 2; k++) {
        y[k] = y0[k] + w[k] * y1[k];
        y[k + n / 2] = y0[k] - w[k] * y1[k];
    }
    return y;
}
  
// Driver code
int main()
{
    vector a{1, 2, 3, 4};
    vector b = fft(a);
    for (int i = 0; i < 4; i++)
        cout << b[i] << endl;
  
    return 0;
}

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

插值
交换 a 和 y 的角色。
用 wn^-1 替换 wn。
将结果的每个元素除以 n。
时间复杂度：O(nlogn)。

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