Macaulay2 » Documentation
Packages » ForeignFunctions » fast Fourier transform example
next | previous | forward | backward | up | index | toc

# 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.

    i2 : libfftw = openSharedLibrary "fftw3"

    o2 = fftw3

    o2 : SharedLibrary

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.

    i3 : fftwPlanDft1d = foreignFunction(libfftw, "fftw_plan_dft_1d", voidstar,
             {int, voidstar, voidstar, int, uint})

    o3 = fftw3::fftw_plan_dft_1d

    o3 : ForeignFunction

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)`.

    i4 : fftwExecute = foreignFunction(libfftw, "fftw_execute", void, voidstar)

    o4 = fftw3::fftw_execute

    o4 : ForeignFunction

    i5 : fftwDestroyPlan = foreignFunction(libfftw, "fftw_destroy_plan", void, voidstar)

    o5 = fftw3::fftw_destroy_plan

    o5 : ForeignFunction

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.

    i6 : fftHelper = (x, sign) -> (
             n := #x;
             inptr := getMemory(2 * n * size double);
             outptr := getMemory(2 * n * size double);
             p := fftwPlanDft1d(n, inptr, outptr, sign, 64);
             registerFinalizer(p, fftwDestroyPlan);
             dbls := splice apply(x, y -> (realPart numeric y, imaginaryPart numeric y));
             apply(2 * n, i -> *(value inptr + i * size double) = double dbls#i);
             fftwExecute p;
             r := ((2 * n) * double) outptr;
             apply(n, i -> value r_(2 * i) + ii * value r_(2 * i + 1)));

    i7 : fastFourierTransform = method();

    i8 : fastFourierTransform List := x -> fftHelper(x, -1);

    i9 : inverseFastFourierTransform = method();

    i10 : inverseFastFourierTransform List := x -> fftHelper(x, 1);

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.

    i11 : fftMultiply = (f, g) -> (
              m := first degree f;
              n := first degree g;
              a := fastFourierTransform splice append(
                  apply(m + 1, i -> coefficient(x^i, f)), n:0);
              b := fastFourierTransform splice append(
                  apply(n + 1, i -> coefficient(x^i, g)), m:0);
              c := inverseFastFourierTransform apply(a, b, times);
              sum(m + n + 1, i -> c#i * x^i)/(m + n + 1));

    i12 : R = CC[x];

    i13 : f = x - 1;

    i14 : g = x^2 + x + 1;

    i15 : fftMultiply(f, g)

           3
    o15 = x  - 1

    o15 : R

The source of this document is in /build/macaulay2-88fgJW/macaulay2-1.25.11+ds/M2/Macaulay2/packages/ForeignFunctions.m2:2198:0.