Algorithms

This page goes from the mathematics of the DFT to what FFTL runs. Each transform in the library is built from a few pieces, and the same pieces appear in the complex, real and real-to-real paths.

Stockham chains

A transform of length   is a sequence of stages. Stage combines interleaved sub-transforms with twiddle factors and writes its result to a second buffer, and the two buffers swap roles at each stage. This is the Stockham formulation, which keeps the data in natural order throughout:

• no bit reversal or digit reversal at the start or the end,

• each stage reads and writes with regular strides,

• the radix sequence is a single object that the plan carries.

The radices are chosen by a fixed rule (src/factor.jl). The rule is the result of exhaustive searches over the decompositions of many lengths on the development machines, not of a runtime search. It prefers radix 4 and radix 8 for powers of two. On aarch64 it merges a factor 3 with a 4 into a 12; the x86 rules differ and are in the same file. Every other prime up to 31 is a radix of its own. Larger primes are not radices; the part of the length they make up goes to Bluestein's algorithm (below).

Codelets

A stage of radix is a small DFT applied to values at once. For each the codelet is straight-line code without loops, generated by a symbolic generator (src/codegen.jl). The generator writes the DFT as a sum of terms and simplifies it:

• multiplications by are rotations of the term, not multiplications; two terms with rotations that differ by an odd power of combine in one complex add-and-rotate;

• a composite radix   is split in registers with the usual index mapping;

• an odd prime uses the symmetric form:    with the cosine sums and the sine sums , which pairs the outputs and  .

The generated code runs on complex SIMD vectors (src/complexvec.jl), with the vector width of the element type. On aarch64 the complex add-and-rotate and the complex multiply are single instructions (FCADD, FCMLA).

Pruned codelets

When the input is known to be zero outside some range, or only some outputs are used, the generator drops the terms that are zero and the statements that feed no needed output. A pass that transforms a zero-padded array therefore does not read the padding, and an NFFT that needs only part of the spectrum does not compute the rest. See transforms along dimensions.

Bluestein's algorithm

Lengths with a prime factor above 31 go through Bluestein's identity. With     ,

The sum is a cyclic convolution with a kernel that depends only on . Its length is a smooth number with   , so the convolution is two transforms of smooth lengths and one pointwise product. The kernel's transform is computed once when the plan is made.

Per line, Bluestein costs two transforms of length and a pointwise product, so it is used only for the part of a length that the smooth codelets cannot cover. On aarch64 the convolution length is the smallest 5-smooth number with   . On x86 a cost model fitted to measured chains picks the fastest 7-smooth length in     : a length with many radix 9, 5 or 3 stages transforms several times slower than a power of two a few percent larger (src/plan.jl).

On x86 a convolution of 8 MiB or more (  in Float64) runs as a four-step transform on the calling thread instead of a chain. A chain of that length no longer fits in cache and makes one pass over memory per stage; the four-step makes two. On a Zen 2 core this is 1.06× faster at and 1.3× from on. On the M1 the chain is faster at every length measured (up to ), so there the chains stay.

Lengths with a smooth part

A length    whose large prime factors make up and whose small ones make up   does not need a convolution of length . Since and are coprime, one Cooley–Tukey step with     and     splits the transform:

1. A chain of length over the interleaved lines (the inner sum),

2. multiplication by the twiddle factors , folded into the first chirp product of the next step,

3. Bluestein transforms of length (the outer sum), each written with stride .

The arithmetic is about the same as one Bluestein transform of length : convolutions of length about are again about points. The gain is where the data lives. A convolution of length fits in cache, one of length often does not. On a Zen 2 core (one thread, Float64)    takes 6.6 ms instead of 16.8 ms, and    takes 43 ms instead of 89 ms.

When itself is large it is split again: becomes the full power of its largest prime, and the first step of length   is again a split transform (or a Bluestein transform when has no smooth part).     runs as Bluestein transforms of length 13213 and 277. Each split costs one more pass over memory, so a leftover stays one Bluestein transform as long as its convolution takes at most 4 MiB. Below that, splitting measured 1.15–1.5× slower; above it, 1.02–1.4× faster (composite_split in src/plan.jl).

Rader's algorithm, which turns a prime length into a cyclic convolution of length  , would halve the convolution where   is smooth. Of 45 random lengths with a large prime factor between and , five qualify; FFTL does not implement it.

Multidimensional transforms without transposition

A transform along dimension of an array with strides is a set of lines with stride . FFTL does not copy the data into lines. The adjacent lines of a dimension with stride are processed together, one SIMD lane per line, so no transposition is needed. Dimensions 1 and 2 are fused per plane when the plane fits in cache.

Long one-dimensional transforms: four-step

A single long line with a smooth length   on several threads uses the four-step algorithm:

1. transforms of length over the columns,

2. multiplication by the twiddle factors ,

3. transforms of length over the rows.

Each of the two transform steps is split over the threads, with a barrier in between. The factorisation is chosen when the plan is made (fourstep_split).

Real transforms

For even  , a real line is read as complex values    , which is the memory layout as it is. One complex DFT of length then gives the half spectrum:

The pairs and   share their products, so the work is half of a complex transform of length . Odd lengths use a complex DFT of length directly. The inverse (c2r) runs the same steps backwards.

Real-to-real transforms

Each kind is computed from a real or complex DFT of a rearranged input. Examples:

• DCT-II. The even-indexed inputs go forward and the odd-indexed ones backward:  . A real DFT of , multiplied by , gives the DCT-II. This is the reordering that makes the DCT cost one real DFT of length .

• DST-II and DST-IV are DCTs of the input with alternating signs, read in reverse order.

• DCT-I and DST-I extend the input to an even or odd sequence of length   or   and take one real DFT of that length.

• DHT is the real DFT with the real and imaginary parts combined:  .

The formulas are in the header comment of src/r2r.jl. Their inverses and scale factors are listed on the cosine, sine and Hartley page.

Passes

plan_pass applies one 1D transform to all lines of a strided array with the same kernel. It lets a library assemble its own multidimensional transform one dimension at a time and declare, per dimension, which inputs are zero and which outputs it needs; both declarations use the pruned codelets above. NFFT.jl transforms its oversampled grid this way.