Fast Fourier Transformation for polynomial multiplication

Given two polynomial A(x) and B(x), find the product C(x) = A(x)*B(x). There is already an O(n^2        ) naive approach to solve this problem. here. This approach uses the coefficient form of the polynomial to calculate the product.
A coefficient representation of a polynomial A(x)=\sum_{j=0}^{n-1}a_jx^j        is a = a0, a1, ..., an-1.

Example-
A(x) = 6x^3 + 7x^2 - 10x + 9        
B(x) = -2x^3 + 4x - 5        
Coefficient representation of A(x) = (9, -10, 7, 6)
Coefficient representation of B(x) = (-5, 4, 0, -2)

Input :
 A[] = {9, -10, 7, 6}
 B[] = {-5, 4, 0, -2}
Output : 
-12x^6 - 14x^5 + 44x^4 - 20x^3 -75x^2 + 86x - 45

We can do better, if we represent the polynomial in another form.

Idea is to represent polynomial in point-value form and then compute the product. A point-value representation of a polynomial A(x) of degree-bound n is a set of n point-value pairs is{ (x0, y0), (x1, y1), ..., (xn-1, yn-1)} such that all of the xi are distinct and yi = A(xi) for i = 0, 1, ..., n-1. 

Example

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

Point-value representation of above polynomial is { (0, 1), (1, 0), (2, 5), (3, 22) }. Using Horner’s method, (discussed here), n-point evaluation takes time O(n^2        ). It's just calculation of values of A(x) at some x for n different points, so time complexity is O(n^2        ). Now that the polynomial is converted into point value, it can be easily calculated C(x) = A(x)*B(x) again using horner's method. This takes O(n) time. An important point here is C(x) has degree bound 2n, then n points will give only n points of C(x), so for that case we need 2n different values of x to calculate 2n different values of y. Now that the product is calculated, the answer can to be converted back into coefficient vector form. To get back to coefficient vector form we use inverse of this evaluation. The inverse of evaluation is called interpolation. Interpolation using Lagrange's formula gives point value-form to coefficient vector form of the polynomial.Lagrange's formula is - 
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})}        
So far we discussed,

.
This idea still solves the problem in O(n^2        ) time complexity. We can use any points we want as evaluation points, but by choosing the evaluation points carefully, we can convert between representations in only O(n log n) time. If we choose “complex roots of unity” as the evaluation points, we can produce a point-value representation by taking the discrete Fourier transform (DFT) of a coefficient vector. We can perform the inverse operation, interpolation, by taking the “inverse DFT” of point-value pairs, yielding a coefficient vector. Fast Fourier Transform (FFT) can perform DFT and inverse DFT in time O(nlogn).

DFT 
DFT is evaluating values of polynomial at n complex nth roots of unity \omega ^{0}_{n},\omega ^{1}_{n},\omega ^{2}_{n}......\omega ^{n-1}_{n}        . So, for y_{k}=\omega ^{k}_{n}        k = 0, 1, 2, ..., n-1, y = (y0, y1, y2, ..., yn-1) is Discrete fourier Transformation (DFT) of given polynomial.
The product of two polynomials of degree-bound n is a polynomial of degree-bound 2n. Before evaluating the input polynomials A and B, therefore, we first double their degree-bounds to 2n by adding n high-order coefficients of 0. Because the vectors have 2n elements, we use “complex 2nth roots of unity, ” which are denoted by the W2n (omega 2n). We assume that n is a power of 2; we can always meet this requirement by adding high-order zero coefficients.

FFT
Here is the Divide-and-conquer strategy to solve this problem.
Define two new polynomials of degree-bound n/2, using even-index and odd-index coefficients of A(x) separately
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)        
The problem of evaluating A(x) at \omega ^{0}_{n},\omega ^{1}_{n},\omega ^{2}_{n}......\omega ^{n-1}_{n}        reduces to evaluating the degree-bound n/2 polynomials A0(x) and A1(x) at the points 
(\omega ^{0}_{n})^{2},(\omega ^{1}_{n})^{2},......(\omega ^{n-1}_{n})^{2}        
Now combining the results by A(x) = A0(x^2) + xA1(x^2) 

The list (\omega ^{0}_{n})^{2},(\omega ^{1}_{n})^{2},......(\omega ^{n-1}_{n})^{2}        does not contain n distinct values, but n/2 complex n/2th roots of unity. Polynomials A0 and A1 are recursively evaluated at the n complex nth roots of unity. Subproblems have exactly the same form as the original problem, but are half the size. So recurrence formed is T(n) = 2T(n/2) + O(n), i.e complexity 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

Pseudo code of recursive FFT

Recursive_FFT(a){
n = length(a) // a is the input coefficient vector
if n = 1
  then return a

// wn is principle complex nth root of unity.
wn = e^(2*pi*i/n)
w = 1

// even indexed coefficients
A0 = (a0, a2, ..., an-2 )

// odd indexed coefficients
A1 = (a1, a3, ..., an-1 ) 

y0 = Recursive_FFT(A0) // local array
y1 = Recursive-FFT(A1) // local array

for k = 0 to n/2 - 1

  // y array stores values of the DFT 
  // of given polynomial. 
  do y[k] = y0[k] + w*y1[k]  
     y[k+(n/2)] = y0[k] - w*y1[k]
     w = w*wn
return y
}
Recursion Tree of Above Execution-

Why does this work?
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}))        
since, \omega ^{k}_{n/2} = \omega ^{2k}_{n} , \omega ^{2k+n}_{n} = \omega ^{2k}_{n} , \omega ^{k+(n/2)}_{n} = -\omega ^{k}_{n}        
Thus, the vector y returned by Recursive-FFT is indeed the DFT of the input
vector a.

#include <bits/stdc++.h>
using namespace std;

// For storing complex values of nth roots
// of unity we use complex<double>
typedef complex<double> cd;

// Recursive function of FFT
vector<cd> fft(vector<cd>& a)
{
    int n = a.size();

    // if input contains just one element
    if (n == 1)
        return vector<cd>(1, a[0]);

    // For storing n complex nth roots of unity
    vector<cd> w(n);
    for (int i = 0; i < n; i++) {
        double alpha = -2 * M_PI * i / n;
        w[i] = cd(cos(alpha), sin(alpha));
    }

    vector<cd> 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<cd> y0 = fft(A0); 

    // Recursive call for odd indexed coefficients
    vector<cd> y1 = fft(A1);

    // for storing values of y0, y1, y2, ..., yn-1.
    vector<cd> 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<cd> a{1, 2, 3, 4};
    vector<cd> b = fft(a);
    for (int i = 0; i < 4; i++)
        cout << b[i] << endl;

    return 0;
}

from math import sin,cos,pi

# Recursive function of FFT
def fft(a):

    n = len(a)

    # if input contains just one element
    if n == 1:
        return [a[0]]

    # For storing n complex nth roots of unity
    theta = -2*pi/n
    w = list( complex(cos(theta*i), sin(theta*i)) for i in range(n) ) 
    
    # Separe coefficients
    Aeven = a[0::2]
    Aodd  = a[1::2]

    # Recursive call for even indexed coefficients
    Yeven = fft(Aeven) 

    # Recursive call for odd indexed coefficients
    Yodd = fft(Aodd)

    # for storing values of y0, y1, y2, ..., yn-1.
    Y = [0]*n
   
    middle = n//2
    for k in range(n//2):
        w_yodd_k  = w[k] * Yodd[k]
        yeven_k   =  Yeven[k]
        
        Y[k]          =  yeven_k  +  w_yodd_k
        Y[k + middle] =  yeven_k  -  w_yodd_k
    
    return Y


# Driver code
if __name__ == '__main__':

    a = [1, 2, 3, 4]
    b = fft(a)
    for B in b:
        print(B)

// JavaScript program to implement
// the approach

// For storing complex values of nth roots
// of unity we use complex
class complex
{
    constructor(a, b = 0)
    {
        this.x = a;
        this.y = b;
    }
    
}

function product(c1, c2)
{
    let c =new complex(0, 0);
    c.x = c1.x * c2.x - c1.y * c2.y
    c.y = c1.x * c2.y + c2.x * c1.y
    return c
}

function sum(c1, c2)
{
    let c = new complex(0, 0);
    c.x = c1.x + c2.x
    c.y = c1.y + c2.y
    return c
}

function difference(c1, c2)
{
    let c =new complex(0, 0);
    c.x = c1.x - c2.x
    c.y = c1.y - c2.y
    return c
}



// Recursive function of FFT
function fft(a)
{
    let n = a.length;

    // if input contains just one element
    if (n == 1)
        return [a[0]]

    // For storing n complex nth roots of unity
    let w = new Array(n);
    let alpha = -2 * Math.PI / n;
    for (var i = 0; i < n; i++) {
        w[i] = new complex(Math.cos(alpha * i), Math.sin(alpha * i));
    }

    let A0 = new Array(Math.floor(n / 2));
    let A1 = new Array(Math.floor(n / 2));
    for (var i = 0; i < Math.floor(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
    let y0 = fft(A0); 

    // Recursive call for odd indexed coefficients
    let y1 = fft(A1);

    // for storing values of y0, y1, y2, ..., yn-1.
    let y = new Array(n);

    for (var k = 0; k < Math.floor(n / 2); k++) {
        
        y[k] =  sum(y0[k], product(w[k], y1[k]));
        y[k + Math.floor(n / 2)] = difference(y0[k], product(w[k], y1[k]));
    }
    return y;
}

// Driver code

let a = [ new complex(1, 0), new complex(2, 0), new complex(3, 0), new complex(4, 0)];
let b = fft(a);
for (var b0 of b)
    console.log("(", b0.x, ",", b0.y, ")")


// This code is contributed by phasing17

import java.util.*;

// For storing complex values of nth roots
// of unity we use complex<double>
class cd {
    double re;
    double im;

    cd(double r, double i) {
        re = r;
        im = i;
    }

    cd add(cd b) {
        return new cd(re + b.re, im + b.im);
    }

    cd sub(cd b) {
        return new cd(re - b.re, im - b.im);
    }

    cd mul(cd b) {
        return new cd(re * b.re - im * b.im, re * b.im + im * b.re);
    }
}

// Recursive function of FFT
class FFT {
    int n;
    cd[] a, w;

    FFT(int n) {
        this.n = n;
        a = new cd[n];
        w = new cd[n];
    }

    void fft() {
        // if input contains just one element
        if (n == 1)
            return;

        for (int i = 0; i < n; i++) {
            double alpha = -2 * Math.PI * i / n;
            w[i] = new cd(Math.cos(alpha), Math.sin(alpha));
        }

        cd[] A0 = new cd[n / 2];
        cd[] A1 = new cd[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
        FFT y0 = new FFT(n / 2);
        y0.a = A0;
        y0.fft();

        // Recursive call for odd indexed coefficients
        FFT y1 = new FFT(n / 2);
        y1.a = A1;
        y1.fft();

        // for storing values of y0, y1, y2, ..., yn-1.
        cd[] y = new cd[n];
        for (int k = 0; k < n / 2; k++) {
            y[k] = y0.a[k].add(w[k].mul(y1.a[k]));
            y[k + n / 2] = y0.a[k].sub(w[k].mul(y1.a[k]));
        }
        System.arraycopy(y, 0, a, 0, n);
    }
// }

// Driver code
// public class Main {
    public static void main(String[] args) {
        FFT fft = new FFT(4);
        fft.a[0] = new cd(1, 0);
        fft.a[1] = new cd(2, 0);
        fft.a[2] = new cd(3, 0);
        fft.a[3] = new cd(4, 0);
        fft.fft();

        for (int i = 0; i < 4; i++)
            System.out.println(fft.a[i].re + " + " + fft.a[i].im + "i");
    }
}

using System;

// For storing complex values of nth roots
// of unity we use a struct named Complex
struct Complex
{
    public double Re;
    public double Im;

    public Complex(double r, double i)
    {
        Re = r;
        Im = i;
    }

    public Complex Add(Complex b)
    {
        return new Complex(Re + b.Re, Im + b.Im);
    }

    public Complex Sub(Complex b)
    {
        return new Complex(Re - b.Re, Im - b.Im);
    }

    public Complex Mul(Complex b)
    {
        return new Complex(Re * b.Re - Im * b.Im,
                           Re * b.Im + Im * b.Re);
    }
}

// Recursive function of FFT
class FFT {
    int n;
    Complex[] a, w;
  
    public FFT(int n)
    {
        this.n = n;
        a = new Complex[n];
        w = new Complex[n];
    }

    public void Fft()
    {
        // if input contains just one element
        if (n == 1)
            return;

        for (int i = 0; i < n; i++) {
            double alpha = -2 * Math.PI * i / n;
            w[i] = new Complex(Math.Cos(alpha),
                               Math.Sin(alpha));
        }

        Complex[] A0 = new Complex[n / 2];
        Complex[] A1 = new Complex[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
        FFT y0 = new FFT(n / 2);
        y0.a = A0;
        y0.Fft();

        // Recursive call for odd indexed coefficients
        FFT y1 = new FFT(n / 2);
        y1.a = A1;
        y1.Fft();

        // for storing values of y0, y1, y2, ..., yn-1.
        Complex[] y = new Complex[n];
        for (int k = 0; k < n / 2; k++) {
            y[k] = y0.a[k].Add(w[k].Mul(y1.a[k]));
            y[k + n / 2] = y0.a[k].Sub(w[k].Mul(y1.a[k]));
        }
        Array.Copy(y, a, n);
    }

    public static void Main()
    {
        FFT fft = new FFT(4);
        fft.a[0] = new Complex(1, 0);
        fft.a[1] = new Complex(2, 0);
        fft.a[2] = new Complex(3, 0);
        fft.a[3] = new Complex(4, 0);
        fft.Fft();

        // css
        // Copy code
        for (int i = 0; i < 4; i++)
            Console.WriteLine(fft.a[i].Re + " + "
                              + fft.a[i].Im + "i");
    }
}

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

Interpolation 
Switch the roles of a and y.
Replace wn by wn^-1.
Divide each element of the result by n.
Time Complexity: O(nlogn).