fast Fourier transform example -- fast Fourier transforms using foreign function calls to the FFTW library

In this example, we make foreign function calls using the FFTW library to compute fast Fourier transforms for polynomial multiplication.

Once FFTW is installed, there should be a file named libfftw3.so or libfftw3.dylib available on the system in one of the standard shared library directories. We can then open the library using openSharedLibrary.

Next, we create our foreign functions using the constructor method foreignFunction. We will need three: one to wrap fftw_plan_dft_1d, which creates the plan for a fast Fourier transform computation, one to wrap fftw_execute, which actually does the computation, and one to wrap fftw_destroy_plan, which deallocates the memory allocated by fftw_plan_dft_1d.

Let's walk through each of these.

According to the FFTW documentation, fftw_plan_dft_1d has the following signature: fftw_plan fftw_plan_dft_1d(int n0, fftw_complex *in, fftw_complex *out, int sign, unsigned flags). Our goal is to compute the discrete Fourier transform $\mathscr F\{\mathbf a\}$ of a vector $\mathbf a\in\CC^n$. In this case, n0 will be $n$, in will be $\mathbf a$ represented as an array of $2n$ double objects (alternating between real and imaginary parts), out will be an array ultimately containing $\mathscr F\{\mathbf a\}$, sign will be 1 unless we want to compute the inverse discrete Fourier transform $\mathscr F^{-1}\{\mathbf a\}$, in which case it will be $-1$. Finally, flags will contain planning flags for the computation. When setting up our foreign function objects, we won't worry about the structures of the fftw_plan or fftw_complex types and will just use voidstar to indicate that they are pointers.

Next, we define foreign functions for fftw_execute, which has signature void fftw_execute(const fftw_plan plan), and fftw_destroy_plan, which has signature void fftw_destroy_plan(fftw_plan plan).

Now we wrap these functions inside Macaulay2 functions. Since computing $\mathscr F\{\mathbf a\}$ and $\mathscr F^{-1}\{\mathbf a\}$ using FFTW requires only changing the sign parameter of the call to fftw_plan_dft_1d, we create a helper function that will do most of the work.

This helper function takes a list x representing the vector $\mathbf a$ as well as the value of sign, allocates the memory required for the in and out lists using getMemory, and then calls the foreign function fftwPlanDft1d. Note that we pass 64 as the value of flags. This specifies the FFTW_ESTIMATE planning flag, which uses a simple heuristic to quickly pick a plan for the computation. We then call registerFinalizer(ForeignObject,Function) to make sure that the fftw_plan object that was returned is properly deallocated during garbage collection. Next, we set up the in array by finding the real and imaginary parts of the elements of x and assigning them using * voidstar = Thing. Finally, we call fftw_execute and convert the result back to a list of complex numbers.

One application of fast Fourier transforms is log-linear time multiplication of polynomials. Indeed, if $f,g\in\CC[x]$ have degrees $m$ and $n$, respectively, and if $\mathbf a,\mathbf b\in\CC^{m+n+1}$ are the vectors formed by taking $a_i$ to be the coefficient of $x^i$ in $f$ and $b_i$ the coefficient of $x_i$ in $g$, then $\frac{1}{m+n+1}\mathscr F^{-1}\left\{\mathscr F\{\mathbf a\}\cdot \mathscr F\{\mathbf b\}\right\}$, where $\mathscr F\{\mathbf a\}\cdot \mathscr F\{\mathbf b\}$ is computed using component-wise multiplication, contains the coefficients of the product $fg$. In some cases, this process can be faster than Macaulay2's native polynomial multiplication.

Let's look at a particular example.