Fast Fourier Transform

Fast Fourier Transform

This is a classic problem in CS – to multiple two degree-n polynomials. The brute force way is obvious, but would cost us O(n²) time. Using the idea of FFT would enable us to do this in O(nlogn) time.
In this blog entry, I will be focusing on the theoretic proof of the FFT algorithm that is so widely used, serving to deepen my own understanding about the topic. Of course, it will be heartening to see if anyone find this article helpful in any way.

Prerequisite knowledge

Polynomials in point value form

For a polynomial of degree-n f(x)=${\sum }_{i=0}^n{a}_i{x}^i$, we choose n arbitrary values of x and use the n pairs of (x_i,f(x_i)) to represent the polynomial. Recall from the gauss’s system of linear equations, the n pair of values is just sufficient for us to figure out what the original polynomial is exactly.

You may have noticed that using the point value form of polynomial greatly reduced the number of operations required for multiplication. Suppose we need to calculate F(x) which is the product of two polynomials f(x) and g(x),then if we have the n pairs of (x_i,f(x_i)) and (x_i,g(x_i)), we can easily derive the corresponding point value form of F(x) which would be (x_i,f(x_i)g(x_i)) in O(n) time. But at what cost? Computing the n values of f(x) alone would require O(n²) time, and we still need F(x) in the coefficient form, and solving the n simultaneous equation would require us O(n³) time by Gauss’s method.

Root of unity

However, notice that we are free to choose any n arbitrary value of x to substitute into the polynomials, and of course we are to choose some special values. In this case, we are to choose n evenly spaced out point on the unit circle of the complex plane, which are the n complex roots of the equation xⁿ=1(recall that in the multiplication of complex numbers can be seen as rotation of the point about the origin on the complex plane, and these n points after rotating n times by angle equal to their argument will end up on the positive side of the real axis with modulus of 1).Let ${\omega }_n=e^{\frac{2\pi i}{n}}$, it is not hard to see this is the first root above the positive real axis, and we can express the whole set of the n roots as
{${\omega }_n^k\mid k=0,1\cdots ,n-1$}, we call them the n^th roots of unity, and they satisfy these beautiful symmetric properties:

1. ${\omega }_n^k={\omega }_n^{k+n}$(going another full circle will not affect anything)
2. ${\omega }_n^k={\omega }_{2n}^{2k}$ (dividing the unit circle into n equal portions and take the k^th crack is equivalent to the 2k^th crack if the unit circle is divided into 2n equivalent portions instead)
3. $-{\omega }_{2n}^k={\omega }_{2n}^{k+n}$(taking the negative value is just reflecting the point through the origin, and the two would differ by an angle of $\pi$, thus we need $\frac{n}{2n}$ of the circle to make same equal)

We are to use these n values in our transformation.

Discrete Fourier Transform

The first step is to convert our f(x) and g(x) from coefficient form into point value form. Take f(x) as an example, the idea is to separate the odd and even powers, and use the idea of Divide and Conquer.

For example, for a polynomial with degree n-1(we assume n is a power of 2 here, as for the missing power we can just assume the coefficient is 0)
$$\begin{array}{rl}\text{suppose we have }f(x)=& \sum _{i=0}^n{a}_i{x}^i\\ \text{let }odd(x)=& \sum _{i=0}^{\frac{n}{2}}{a}_{2\ast i+1}{x}^i,\\ even(x)=& \sum _{i=0}^{\frac{n}{2}}{a}_{2\ast i}{x}^i,\\ \text{Then we have }f(x)=& even(x^2)+x\ast odd(x^2)\end{array}$$
Now we can substitute in our n^th roots of unity:
$$\begin{array}{rl}f({\omega }_n^k)& =even(({\omega }_n^k{)}^2)+{\omega }_n^k\ast odd(({\omega }_n^k{)}^2)\\ & =even({\omega }_n^{2k})+{\omega }_n^k\ast odd({\omega }_n^{2k})\\ & =even({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\ast odd({\omega }_{\frac{n}{2}}^k)\text{...from(2)}\\ f({\omega }_n^{k+\frac{n}{2}})& =even(({\omega }_n^{k+\frac{n}{2}}{)}^2)+{\omega }_n^{k+\frac{n}{2}}\ast odd(({\omega }_n^{k+\frac{n}{2}}{)}^2)\\ & =even({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}\ast odd({\omega }_n^{2k+n})\\ & =even({\omega }_n^{2k})-{\omega }_n^k\ast odd({\omega }_n^{2k})\text{...from(1),(3)}\\ & =even({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\ast odd({\omega }_{\frac{n}{2}}^k)\text{...from(2)}\end{array}$$
Thus, we only need to compute $even({\omega }_{\frac{n}{2}}^k)$ and $odd({\omega }_{\frac{n}{2}}^k)$ to derive the values of $f({\omega }_n^k)$ and $f({\omega }_n^{k+\frac{n}{2}})$, and we can carry out this process recursively for the odd and even functions as well. For each recursion, we convert the computation of a degree-n polynomial into 2 degree-$\frac{n}{2}$ polynomials, and we only compute half the number of point values (${\omega }_n^k$) for both the 2 degree-$\frac{n}{2}$ polynomials since the other half (${\omega }_n^{k+\frac{n}{2}}$) just differ by a plus/minus sign. Thus, we need $2\ast \frac{1}{2}\ast n=n$ operations for each round of recursion, and after $lo{g}_{2}n$ rounds, we will end up with n degree-1 polynomials and we only need to substitute 1 into each one of them. This would give us a time complexity of O$(nlogn)$.

This process is known as the Discrete Fourier Transform(DFT).

Inverse Discrete Fourier Transform

The second part of the FFT algorithm would be reverting back from point value form to coefficient form. As you might have guessed, the process of Inverse Discrete Fourier Transform(IDFT) is largely similar to that of DFT.

In fact, if we see DFT as a linear transformation from the vector of coefficients {$a_0,a_1,a_2,...,a_{n-1}$} to the vector of the n point values {$f(x_0),f(x_1),f(x_2),...,f(x_{n-1})$}, we are essentially multiplying the coefficient vector by an n$\times$n matrix. In fact, the brute force way of substituting each of the n values one by one is essentially matrix multiplication, and our DFT just simplifies the process by using some special properties of root of unity(and by Divide and Conquer).
$$\left[\begin{array}{c}f(x_0)\\ f(x_1)\\ f(x_2)\\ f(x_3)\\ ⋮\\ f(x_{n-1})\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \cdots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \cdots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]$$
Conversely, if we have the point value vector, and we want to compute the coefficient vector, we just need to multiply the inverse matrix of the bulky $n\times n$ matrix on the left for both sides.
Again, inverse matrices are generally slow to compute, but in this case we can check that the following matrix is the inverse matrix.
$$\frac{1}{n}\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& {\omega }_n^{-1}& {\omega }_n^{-2}& {\omega }_n^{-3}& \cdots & {\omega }_n^{-(n-1)}\\ 1& {\omega }_n^{-2}& {\omega }_n^{-4}& {\omega }_n^{-6}& \cdots & {\omega }_n^{-2(n-1)}\\ 1& {\omega }_n^{-3}& {\omega }_n^{-6}& {\omega }_n^{-9}& \cdots & {\omega }_n^{-3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{-(n-1)}& {\omega }_n^{-2(n-1)}& {\omega }_n^{-3(n-1)}& \cdots & {\omega }_n^{-(n-1{)}^2}\end{array}\right]$$
We can easily check that the product of the two matrices is the identity matrix. Thus, the IDFT process would be very similar to DFT except for we are using the “negative” roots of unity now. In fact, we can reuse the code of DFT for IDFT.

Sample Code

#include <bits/stdc++.h>
using namespace std;
#define LL long long
#define il inline
const int N=3e5+5;
const double pi=acos(-1);
struct comp//complex number
{
    double x,y;
    comp operator + (const comp &c)const{
        return {x+c.x,y+c.y};
        }
    comp operator - (const comp &c)const{
        return {x-c.x,y-c.y};
    }
    comp operator * (const comp &c)const{
        return {x*c.x-y*c.y,x*c.y+y*c.x};
    }
}a[N],b[N];
int n,m,tot=1,bit=0,rev[N];
void fft(comp c[],int inv)//IFFT if inv=1
{
    for(int i=0;i<tot;i++){
        if(i<rev[i])swap(c[i],c[rev[i]]);
    }
    for(int mid=1;mid<tot;mid<<=1){
        comp w1={cos(pi/mid),inv*sin(pi/mid)};
        for(int i=0;i<tot;i+=mid*2){
            comp wk={1,0};
            for(int j=0;j<mid;j++,wk=wk*w1){
                comp x=c[i+j],y=wk*c[i+j+mid];
                c[i+j]=x+y,c[i+j+mid]=x-y;
            }
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lf",&a[i].x);
    for(int i=0;i<=m;i++)scanf("%lf",&b[i].x);
    while(tot<n+m+1){bit++,tot<<=1;}
    for(int i=0;i<tot;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    fft(a,1),fft(b,1);
    for(int i=0;i<tot;i++)a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/tot+0.2))//round off;
}