Why is FFTW written in OCaml and what makes it so fast?

https://en.wikipedia.org/wiki/Talk%3AFFTW#Bernstein_%22contr...

I carefully studied the source code of FFTW a few years ago and read a few of Frigo's papers (e.g., [2]). The purpose of FFTW is to apply direct/inverse Fourier transforms to vectors of N elements, and it does so by employing a number of algorithms to split the data into many smaller chunks and then applying the transform to each of these sub-chunks. (This is basically the idea of the FFT.)

When a chunk is small enough (e.g., 10 elements), it can be worth to stop splitting it into smaller subpieces but instead to apply a direct brute-force formula. This is where OCaml comes into the game. FFTW's authors have written OCaml snippets to compute the Fourier Transform in a number of ways. This code is written using custom operators, so that the snippet is not compiled down to machine code but instead kept in a AST. An optimizing complier (written in OCaml) is then ran on the AST to perform a number of optimizations to the snippet, and finally output C code is generated out of the optimized AST. Among the optimizations there are some classical ones (e.g., constant folding); others are domain-specific, because they employ some property or symmetry of the Fourier Transform. The C code that is produced as output is taylored for some specific situation, e.g., direct transform of 13 purely real numbers.

The algorithm that takes N numbers and splits them in chunks is written in C, as well as the code that decides when to stop the recursive splitting and which OCaml-generated snippet is best suited for each chunk. Moreover, C is used for all the function that allocate and handle the many data structures used by the library. Therefore, it is a bit of a stretch to say that FFTW is written in OCaml, although its most interesting part (the optimizing compiler) is indeed. This is the reason why Frigo's paper [2] deals mainly with the OCaml part.

[1] http://www.ffconsultancy.com/products/ocaml_for_scientists/

[2] http://www.fftw.org/fftw-paper-ieee.pdf

FWIW, proper hygienic macros and scientific focused array types makes this sort of meta-programming optimization in Julia kinda fun.

Though just using type dispatching / multi-methods might be all that’s needed. It all makes me want to have a reason to write that kind of optimized code... For example if you can lift the size, variable, and variability into types you could match on your config roughly like: “fft(X :: Pwr2_64, Y :: FWD, Dir :: REV, Var :: BIDIR, Z) = 2^4-2^(log2(X))”.

> FFTW_MEASURE tells FFTW to find an optimized plan by actually computing several FFTs and measuring their execution time. Depending on your machine, this can take some time (often a few seconds). FFTW_MEASURE is the default planning option.

> FFTW_PATIENT is like FFTW_MEASURE, but considers a wider range of algorithms and often produces a “more optimal” plan (especially for large transforms), but at the expense of several times longer planning time (especially for large transforms).

> FFTW_EXHAUSTIVE is like FFTW_PATIENT, but considers an even wider range of algorithms, including many that we think are unlikely to be fast, to produce the most optimal plan but with a substantially increased planning time.

http://www.fftw.org/fftw3_doc/Planner-Flags.html

I believe this sort of technique is used for some cryptography libraries as well, but I don't have any specifics.

The ATLAS BLAS/LAPACK library also uses a similar technique.

i was like

DUDE YOU MUST SEE THIS IS THE COOLEST APPLICATION OF OCAML COMPILERS EVER