Discrete Fourier Transform

Suppose we have to multiply two degree polynomials together, $A(x)B(x) = C(x)$ where $A$ and $B$ are given, $C$ is unknown. Let $A,B,C$ be of degree $r,s,r+s$ respectively, and let $n$ be the smallest power of 2 larger than $r+s$.  \[    A(x) = a_0 + a_1x + a_2x^2 + ... + a_{n-1}x^{n-1}     \] \[    B(x) = b_0 + b_1x + b_2x^2 + ... + b_{n-1}x^{n-1}     \] \[    C(x) = c_0 + c_1x + c_2x^2 + ... + c_{2n-2}x^{2n-2}     \]
Note that the coefficients indexed larger than the respective degrees are simply zero. Then we can find our coefficients of $C$ in $\mathcal{O}(n^2)$ by multiplying the coefficients of $A$ and $B$ together in the convoluted sum \[    c_k = \sum_{j=0}^{k}a_j b_{k-j}    \]
But $n=2^{20}$ and $\mathcal{O}(n^2)$ is too costly, so we need a better method. Through the use of the FFT we can find coefficients $c_k$ in $\mathcal{O}(n\log n)$ time by

1.     (FFT) Evaluate $A$ and $B$ at $n$th roots of unity $\omega^k, \;\;\;\; k=0,...,n-1$
2.      Find $C(\omega^k) = A(\omega^k) B(\omega^k)$  
3.      (iFFT) Find $C$ via polynomial interpolation  

Step (1) Suppose we write $A(x) = A_0(x^2) + xA_1(x^2)$ where \[   A_0(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{n/2-1} \] \[    A_1(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{n/2-1}  \]
To evaluate $A$ at the $n$th roots of unity, we need to evaluate $A_0$ and $A_1$ at the square of each of the roots. But this is the same as evaluating the $(n/2)$th roots of unity at two polynomials of degree $n/2-1$. We immediately get a recursive algorithm (Fast Fourier Transform) with complexity $T(n) = 2T(n/2) + 2n = \mathcal{O}(n\log n).$ We can easily generalize this to work for $n$ as a power of any prime by splitting $A$ into more parts (or any number provided you have its prime factorization).

Step (2) We evaluate $C(\omega^k )$ in $O(n)$ operations.

Step (3) We have formed a linear system $Mc^* = c$, where $M$ is the matrix $(\omega^{kj})$, $c^*$ is the vector of values $ c_k $, and $c$ is the vector of values $C(\omega^{km})$. $M$ is a Vandermonde matrix, so it has an inverse. Thus, we can find $c^{*}$ by applying $M^{-1}$ to both sides, $c^* = M^{-1} c$. It turns out $M^{-1} = (1/n (\omega^{-ij}))$ and applying this matrix to $c$ is the same as step 1, namely the inverse Fourier transform.

Python Code




    from cmath import exp, pi 

    def FFT(X):
       
        n = len(X)
        w = exp(-2*pi*1j/n)

        if n > 1
            X = FFT(X[::2]) + FFT(X[1::2])

            for k in xrange(n/2):

                xk = X[k]
                X[k] = xk + w**k*X[k+n/2]
                X[k+n/2] = xk - w**k*X[k+n/2]

        return X


___________________________________________________________________________________