📜  快速傅立叶变换用于多项式乘法

📅  最后修改于: 2021-05-07 10:02:53             🧑  作者: Mango

给定两个多项式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都是唯一的,并且对于i = 0、1,…,n-1,yi = A(xi)。
例子

 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 )。现在将多项式转换为点值,可以再次使用霍纳法轻松地计算出C(x)= A(x)* B(x)。这需要O(n)时间。这里的一个重要点是C(x)的度界为2n,那么n个点将仅给出n个C(x)点,因此对于这种情况,我们需要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
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个元素,因此我们使用“复数第2个单位的第n个根”,表示为由W2n(Ω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处评估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个复数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)//如果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 =递归FFT(A1)// k = 0到n / 2-1的本地数组// y数组存储DFT //的值给定多项式。 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}
    因此,递归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)。