FFT

Numerical

Produce two polynomials, and their product.

Unprotect[C];
Clear[A, B, C]
A[x_]  :=  (x^3)  +  2 (x^2)  -  x  +  5
B[x_]  :=  2 (x^3)  +  5 (x^2)  -  3 x  -  1
C[x_]  :=  Evaluate@Expand[ A[x] B[x] ]
C[x]

A primitive 8th root of unity in ℂ.

ω  =  E^( 2 I Pi/8)

Evaluating the polynomials.

AA  =  
      Map[ A,  (ω^Range)[0, 7] ]  //  
        Simplify
BB  =  
      Map[ B,  (ω^Range)[0, 7] ]   //  Simplify

Pointwise product.

V1  =  AA BB   //  Simplify

V2 =  Map[ C ,  (ω^Range)[0, 7]]  //  Simplify

Compare the value lists.

V1  -  V2

V1  -  V2  //  Simplify

Symbolic

Produce two symbolic polynomials, and their product.

Clear[A, B, C]
A[x_]  :=  
      (x^3)  +  (Subscript[a, 2]) (x^2)  +  (Subscript[a, 1]) x  +  (Subscript[a, 0])
B[x_]  :=  (x^3)  +  (Subscript[b, 2]) (x^2)  +  (Subscript[b, 1]) x  +  (Subscript[b, 0])
C[x_]  :=  Evaluate@Expand[ A[x] B[x] ]
C[x]

A primitive 8th root of unity in ℂ.

ω  =  E^( 2 I Pi/8)

Evaluating the polynomials.

AA  =  
      Map[ A,  (ω^Range)[0, 7] ]  //  
        Simplify
BB  =  
      Map[ B,  (ω^Range)[0, 7] ]   //  Simplify

Pointwise product.

V1  =  AA BB   //  Simplify

V2 =  Map[ C ,  (ω^Range)[0, 7]]  //  Simplify

Compare the value lists.

V1  -  V2

V1  -  V2  //  Simplify

Note:  if you trust the arithmetic in Mathematica, this constitutes a proof that our trick works for all (monic) polynomials of degree 3.

Recursive FFT

A direct implementation of FFT.

Clear[FFT]

FFT[ {a_} ] := {a};

FFT[ A_List ] := 
Module[ {n, om,  Om, Yeven, Yodd, Aeven, Aodd},

	n  =  Length[A];
	om  =  E^(2 I Pi/n);
	Om  =  om^Range[0, n/2 - 1];

  
	{Aeven,  Aodd}  =  Thread@Partition[ A,  2 ];
	Yeven = FFT[ Aeven ];
	Yodd  = FFT[ Aodd ];

	Join[  Yeven  +  Om * Yodd,  
                Yeven  -  Om * Yodd ]
];

Comparison to system function.

FFT[ Range[0, 7] ]  //  N
Fourier[ Range[0, 7],  
      FourierParameters -> {1, 1} ]

FFT[ Range[2^10] ];  //  Timing
Fourier[ Range[2^10] ];  //  Timing

First[Timing[FFT[ Range[2^#] ]]] &   /@  Range[10]  //  
      TableForm

Clear[f, g, h]
f[x_]  :=  
      x^3  +  2 x^2  -  x  +  5
g[x_]  :=  
      2 x^3  +  5 x^2  -  3 x  -  1
h[x_]  :=  Evaluate@Expand[ f[x] g[x] ]

ff  =  
      PadRight[CoefficientList[f[x], x], 8]
gg  =  
      PadRight[CoefficientList[g[x], x], 8]
hh  =  PadRight[CoefficientList[h[x], x], 8]

Fourier[ ff,  
      FourierParameters -> {1, 1} ]
Fourier[ gg,  
      FourierParameters -> {1, 1} ]
L  =  % * %%

InverseFourier[ L,  FourierParameters -> {1, 1} ]

Inverse FFT

Vandermonde[ a_,  
        n_ ]  :=  
Table[  a^(i j),  {i, 0, 
          n - 1},  {j, 0, n - 1} ]

Vandermonde[x,  6 ] //  MatrixForm

M  =  
        Vandermonde[ E^(2 Pi I/4),  
          4 ];
M  //  MatrixForm

Inverse[M]  //  Simplify  //   MatrixForm

A small example:

Y  =  
      Fourier[ Range[0, 7],  FourierParameters -> {1, 1} ]

The system function:

InverseFourier[ Y,  FourierParameters -> {1, 1} ]

And a simple implementation by hand:

Clear[iFFT]
iFFT[ {a_} ] := {a};
iFFT[ A_List ] := 
Module[ {n, om,  Om, Yeven, Yodd, Aeven, Aodd},
	
	n  =  Length[A];
	om  =  E^(-2 I Pi/n);
	Om  =  om^Range[0, n/2 - 1];

  
	{Aeven,  Aodd}  =  Thread@Partition[ A,  2 ];
	Yeven = iFFT[ Aeven ];
	Yodd   = iFFT[ Aodd ];

	Join[  Yeven  +  Om * Yodd,  
                Yeven  -  Om * Yodd ]
];

iFFT[ Y ]