One chain form: every smooth length runs on RChain

docs/AGENTS.md

@@ -35,8 +35,8 @@ rules there apply to everything written here as well. Everything in this folder
 - **Explanation** (`src/explanation/`): the DFT and its conventions, FFTW and FFTL (history, papers, what
   FFTL takes from FFTW; written from papers and documentation only), algorithms (Stockham chains,
   codelets, Bluestein, multidimensional without transposition, four-step, real and real-to-real
-  reductions), implementation (type specialisation, SIMD, kernel cache, cache sizes, the two chain forms,
-  why libfftl needs `RChain`), design decisions (what was decided and why, with the evidence), package
+  reductions), implementation (type specialisation, SIMD, kernel cache, cache sizes, the stage loop of the
+  chain, why libfftl needs `RChain`), design decisions (what was decided and why, with the evidence), package
   size and load time, performance (measurement protocol).
 - **Reference** (`src/reference/`): Julia API (`@docs`), AbstractFFTs (what is implemented), libfftl C API
   (prose on FFTW3's interface as libfftl implements it; the header pages are generated from
@@ -46,7 +46,7 @@ rules there apply to everything written here as well. Everything in this folder
 
 ## Decisions for the manual
 
-The design decisions of the package (interface, chain forms, rules for lengths, package size) are on the
+The design decisions of the package (interface, the chain, rules for lengths, package size) are on the
 Design decisions page. These concern the manual itself:
 
 - **Interface in the text:** the manual writes the exported names without prefix; next to `using FFTL` it

docs/src/explanation/design-decisions.md

@@ -74,29 +74,49 @@ with one thread, and from ×0.77 to ×1.16 with 16 threads.
 measured on Zen 2 only, and the 4 MiB rule is checked on 12 lengths that also took part in its
 derivation.
 
-## Two chain forms: compiled for short smooth lengths, run-time data for the rest
-
-**Decision (2026-11-09).** `kernel` builds the compiled `Chain` (radix sequence in the type) for smooth
-lengths n ≤ 128 (`SMALL_N`) and the `RChain` (radix sequence as data, a fixed set of codelets) for
-every other smooth length. The choice is made when the plan is built, so a call carries no dispatch.
-Both forms share the twiddle tables, the radix rule and the vector-type rules (bitwise equal in 552
-checks). `RChain` is also the only form the trimmed C library can contain.
-
-**Why.** A plan of a new size with `Chain` is compilation, not planning: 0.31–0.48 s for plan and first
-call against 0.013–0.044 s with the run-time plan (M1, one fresh process per case). But on x86 the
-run-time choice of the codelet costs 9–20 ns per call (n = 16, `Float64`: `Chain` 22 ns, `RChain`
-39 ns), which is 1.2–1.5× for short lines. With the split, short lengths are at parity or better. Warm execution of the run-time plan against the
-compiled one was 1.00–1.08 on the cases measured, with one exception (Bluestein, 100003, `Float32`: 1.42,
-cause open).
-
-**Replaced.** The idea of removing `Chain` altogether (2026-11-08) and the earlier decision to keep the
-JIT and run-time plan families as separate types (2026-11-02).
+## One chain form: the radix sequence is data
+
+**Decision (2026-11-13).** Every smooth length runs on `RChain`: the radix sequence and the stage
+records are data, the codelets come from a fixed set. There is no chain type with the radix sequence
+in its type, and nothing is compiled per length. `RChain` is also the only form the trimmed C library
+can contain.
+
+**Why.** The compiled chain (`Chain{T,Rs,σ,K}`, used for smooth lengths up to 128 until this decision)
+existed because the run-time chain cost 9–20 ns more per call on x86. That cost was not the choice of
+the codelet at run time (about 1 ns per stage) but the code around it: GC frames for references loaded
+from the chain, the copy of an immutable chain into each caller, bounds checks, and one more call level
+per stage. With those removed (see [Implementation](@ref "Implementation")) the compiled chain has no
+advantage left. AMD EPYC 9354 (AVX-512), one thread, `plan_fft` and `mul!`, 92 cases (`Float64` and
+`Float32`; lines of 2 to 4096, planes of 8 × 8 to 128 × 128, cubes of 8³ to 64³, in-place, backward and
+real transforms), 8 rounds in alternating processes, median time over the median time with the compiled
+chain:
+
+- geometric mean 0.973; 27 cases below 0.95, 8 above 1.05;
+- lines of one stage (2, 4, 16): 0.80–0.95; 8 × 8, 8 × 8 × 8, 7 × 12, 5 × 64: 0.85–0.91;
+- the cases above 1.05 are short lines of two or three stages: `Float32` 20, 32, 96 and 100 (1.09–1.10),
+  37 rows of 64 (1.13), `rfft` of 64 (1.07 and 1.09), and `Float64` 12 × 7 (1.05). The stage loop runs
+  about 50 instructions per stage around the codelet (no mispredicted branches), 2.4 ns per stage
+  against direct calls of the two kernels of `Float32` 32.
+- The machine was shared during the measurement: over the 8 rounds the slowest time of the compiled
+  chain over its fastest is 1.04 in the median over the cases, 1.19 at the 90th percentile and 1.55 at
+  most.
+
+Not verified on aarch64.
+
+What the single form removes: 583 compiled chain methods and their generators, the dispatch trees from
+an `RChain` of a short length into the compiled chains, and the compilation at the first plan of a
+short length with a radix sequence outside the compiled set (a short dimension with a large stride).
+The package has 1796 methods instead of 2385, the package image on this machine is 61 MB instead
+of 72 MB, and `libfftl.so` (x86-64-v3, `--trim=safe`) is 49.5 MB instead of 55.6 MB.
+
+**Replaced.** The two chain forms (2026-11-09: compiled `Chain` for smooth lengths up to 128, `RChain`
+for the rest), and before that the idea of removing `Chain` without changing the stage loop
+(2026-11-08) and the separate JIT and run-time plan families (2026-11-02).
 
 ## Kernels of a fixed set are plain methods, not `@generated`
 
-**Decision (2026-11-10).** The stage kernels of the codelet set that `RChain` dispatches to, and the
-compiled chains of the short lengths, are plain methods defined in `src/kernel_methods.jl` from the
-bodies the `@generated` methods build. The `@generated` methods remain as the fallback.
+**Decision (2026-11-10).** The stage kernels of the codelet set that `RChain` dispatches to are plain
+methods defined in `src/kernel_methods.jl` from the bodies the `@generated` methods build. The `@generated` methods remain as the fallback.
 
 **Why.** A package image keeps the complete optimised IR of every precompiled `@generated`
 specialisation, 130 MB of the 146 MB image. See [Package size and load time](@ref "Package size and load time").
@@ -157,9 +177,8 @@ and matches the direct DFT to 3·10⁻¹⁴ when called from C. Planning on 64 
 `d = 1` and 0.9 µs for `d = 2` on the M1 (the specialised plan: 4.7 µs and 2.8 µs). Applying the
 specialised and the run-time plans is within noise of each other on the same views.
 
-**Costs.** The run-time kernels dispatch per stage, which costs a little on short lines (the reason the
-specialised plans stay the default). For real lengths with large prime factors, the run-time kernel is a
-Bluestein transform, without the composite split of the specialised plan.
+**Costs.** For real lengths with large prime factors, the run-time kernel is a Bluestein transform,
+without the composite split of the specialised plan.
 
 ## The interface: FFTW.jl's names, FFTL's own functions, AbstractFFTs as a backend
 

docs/src/explanation/image-size.md

@@ -43,10 +43,9 @@ later Julia versions, so re-run the count below after an upgrade.
 ## What keeps it small
 
 1. **Plain methods for the fixed set.** The stage kernels that `RChain` dispatches to (radix × sign ×
-   vector type, empty zero-input mask or full output mask) and the compiled chains of the smooth lengths
-   up to `SMALL_N` are defined in `src/kernel_methods.jl` with `@eval`, from the same body-building
+   vector type) are defined in `src/kernel_methods.jl` with `@eval`, from the same body-building
    functions the generators call. The `@generated` methods stay as the fallback for every other combination
-   (other radices, vector types and chains), so behaviour is unchanged. A plain method keeps its (compressed) source
+   (the lazy radices), so behaviour is unchanged. A plain method keeps its (compressed) source
    and its machine code, not its optimised IR.
 2. **Lazy radices.** The codelets of the radices in `LAZY_RADICES` (17, 19, 23, 29, 31 by default; the
    preference `lazy_radices`, read in `src/rchain.jl`) are 54 % of the statements of the stage kernels and are rarely needed. They are not
@@ -69,6 +68,10 @@ later Julia versions, so re-run the count below after an upgrade.
 To add a radix, put it in `RADICES`; `kernel_methods.jl` defines its plain methods. To make a radix lazy,
 add it to the default of the preference in `src/rchain.jl`.
 
+The numbers on this page were measured while the smooth lengths up to 128 had compiled chains of their
+own (583 more plain methods). Without them the image on an AMD EPYC 9354 (AVX-512, which has more vector
+types than the M1) went from 72 MB to 61 MB; the M1 numbers have not been measured again.
+
 ## What is left
 
 The 32 MB consist of about 3.8 MB for the package without the kernels, 8.6 MB for the definitions of the

docs/src/explanation/implementation.md

@@ -3,18 +3,16 @@
 The [algorithms](@ref "Algorithms") explain what is computed. This page explains how the code
 is organised and why, for someone who wants to change it.
 
-## Specialisation by type, not by configuration
+## Specialisation by type for the codelets, data for the plan
 
-A plan is a Julia value, and its type says what it is. A Stockham chain is
-`Chain{T,Rs,σ,K}`: the element type `T`, the radix sequence `Rs` as a tuple of integers, the sign
-`σ` and the kernel type `K` are type parameters. The compiler therefore sees the complete stage
-sequence when it compiles the transform, unrolls the loop over stages where it can, and emits
-machine code for exactly that sequence. The price is one compilation per radix sequence, which
-is what the [plan cache](@ref "Plans and reuse") and the [precompile workload](@ref "Precompilation")
-are for.
+A codelet is specialised by type: the stage kernels take the radix, the sign and the vector type as
+type parameters, and the straight-line code of each combination is compiled once. The combinations
+form a fixed set (radix × sign × vector type × element type), so all of them are compiled with the
+package.
 
-The codelets follow the same rule. `@generated` functions build the straight-line code for a
-given radix, element type and vector width when they are first called with those types.
+The stage sequence of a length is data. A Stockham chain (`RChain`, `src/rchain.jl`) holds its
+radix sequence and one record per stage (codelet index, stride, twiddle table), and a loop calls the
+codelet of each stage. Nothing is compiled when a plan of a new length is made.
 
 ## Codelets on SIMD vectors
 
@@ -56,25 +54,28 @@ buffers per task are the minimum: the Stockham stages ping-pong between them.
 That buffer ownership is why a plan must not run on two tasks at once. The plan is a
 container for the buffers, not a stateless function.
 
-## Two forms of the same chain
+## The stage loop
 
-There are two representations of a Stockham chain:
+`RChain` calls codelets from the fixed set in a plain loop, selected by an integer index. A short line
+is a few nanoseconds of arithmetic, so the loop is written to cost little next to it:
 
-- `Chain` (`src/plan.jl`) has the radix sequence in its type and the stage sequence is unrolled at compile
-  time. It is the fastest form for short lines: a call carries no dispatch.
-- `RChain` (`src/rchain.jl`) has the radix sequence as data and calls codelets from a fixed set (radix ×
-  sign × vector type × element type) in a plain loop, selected by an integer index.
+- A line that is one stage is one codelet call, without the loop.
+- The stage records lie in one flat vector and carry the address of their twiddle table, and the loop
+  reads them through a raw pointer. A reference loaded from the chain inside the loop would need a GC
+  frame, which costs more than a short stage.
+- The chain is a mutable struct, so it is passed as one pointer; an immutable chain is copied to the
+  stack of each caller.
+- The switch over the codelet index is inlined into the loop: a stage is one jump and one call.
+- The rows of a plane run inside one call of the loop function.
 
-The plan chooses when it is built: `Chain` for smooth lengths up to `SMALL_N` = 128, `RChain` for every
-other smooth length. A new radix sequence of `Chain` is a new type, so the first plan of such a size is
-compilation; `RChain` has no per-size compilation, and its codelet set is finite, so it can be
-precompiled completely. On x86 the run-time choice of the codelet costs 9–20 ns per call, which is why
-the short lengths keep the compiled form.
+The addresses in the stage records are valid for the chain object they were taken for. A chain
+remembers its own address at that moment; a copy of a plan, or a plan stored in a package image and
+loaded into another process, is at another address and links its records again before it uses them.
 
 `RChain` is also what the C library is built from. `juliac --trim` removes everything it cannot resolve
-at build time, and a chain whose type is only known at run time cannot be compiled that way. The
+at build time and compiles nothing at run time; the fixed codelet set is all a chain of any length needs. The
 [libfftl how-to](@ref "C library (libfftl)") describes what the library is built from, and
-[Design decisions](@ref "Design decisions") gives the measurements behind the split. The stage kernels of
+[Design decisions](@ref "Design decisions") gives the measurements behind the single chain form. The stage kernels of
 the fixed codelet set are plain methods, not `@generated` ones, to keep the package image small; see
 [Package size and load time](@ref "Package size and load time").
 

src/interface.jl

@@ -116,7 +116,7 @@ check_size(p, x) = size(x) == p.sz ||
 # Dimensions 1 and 2 are fused per plane when a plane fits in cache (`run_planes!`).
 function apply_dims!(y::StridedArray{Complex{T}}, x::StridedArray{Complex{T}}, p::FFTLPlan{T}, ks::Tuple, ds::Tuple) where {T}
   isempty(ks) && return nothing
-  if length(ks) >= 2 && ds[1] == 1 && ds[2] == 2 && ks[1] isa AnyChain && ks[2] isa AnyChain
+  if length(ks) >= 2 && ds[1] == 1 && ds[2] == 2 && ks[1] isa RChain && ks[2] isa RChain
     n1, n2 = p.sz[1], p.sz[2]
     outer = dim_extent(p.sz, 2)[2]
     if fuse_planes(n1, n2, outer, p.nthreads)
@@ -193,7 +193,7 @@ function run_strided_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Kernel{T}, d::Int,
   groups = CartesianIndices(ntuple(j -> j <= m || j == d ? 1 : p.sz[j], N))
   # LinesMT owns a kernel and buffers per task; a composite kernel holds at most k.B lanes.
   nt = k isa LinesMT ? min(p.nthreads, length(k.kernels)) :
-       k isa AnyChain && threaded(lanes, n, length(groups), p.nthreads) &&
+       k isa RChain && threaded(lanes, n, length(groups), p.nthreads) &&
        length(p.tbuf[1]) >= n ? p.nthreads : 1
   nw = nt == 1 ? length(p.work1) : length(p.tbuf[1])
   B = k isa LinesMT ? min(lanes, k.B) : min(lanes, lane_block(lanes, n, nt), nw ÷ n)

src/kernel_methods.jl

@@ -1,5 +1,5 @@
 # Compiled-ahead methods of the stage kernels for the fixed codelet set (radix × sign × vector
-# type) that `RChain` dispatches to, and of the chains of the small lengths (rchain.jl). The
+# type) that `RChain` dispatches to. The
 # `@generated` kernels in kernels.jl define every other combination. A package image keeps the
 # complete optimised IR of each precompiled `@generated` specialization but none of a plain
 # method, so the fixed set is defined as plain methods with the same bodies. The codelets of
@@ -48,15 +48,3 @@ for T in (Float32, Float64), σ in (-1, 1), R in RADICES
   b = kernel_subst(stage_small_body(T, R, σ, VT), Dict{Symbol,Any}(:T => T, :VT => VT))
   @eval $(Expr(:function, h, b))
 end
-
-# The compiled chains of the small lengths (rchain.jl) run through `run!`/`run_lanes!` methods of
-# plain definition for the same reason; any other chain type uses the generic methods.
-for T in (Float32, Float64), (n, σ) in SMALL_ENTRIES
-  Rs = Tuple(radices(n, 1; esize = sizeof(Complex{T})))
-  C = Chain{T,Rs,σ,length(Rs)}
-  @eval $(Expr(:function, :(run!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::$C, S::Int)),
-               run_body(T, Rs, σ, length(Rs))))
-  @eval $(Expr(:function, :(run_lanes!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::$C,
-                                       Sx::Int, Sy::Int, B::Int)),
-               run_lanes_body(T, Rs, σ, length(Rs))))
-end

src/kernels.jl

@@ -46,23 +46,6 @@ end
   stage_body(T, R, σ, VT)
 end
 
-# Stage for strides that are not a multiple of the complex lanes of VT (odd batch sizes): the
-# first s - s % C lanes in vectors of C complex values, the scalar remainder one complex value
-# at a time.
-@inline function stage_split!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, s::Int,
-                              ::Val{R}, ::Val{σ}, ::Type{VT}) where {T,R,σ,VT}
-  C = ncomplex(VT)
-  rem = s % C
-  if rem == 0
-    stage!(py, px, pt, m, s, Val(R), Val(σ), VT)
-  else
-    stage_range!(py, px, pt, m, s, 0, s - rem, 0, m, Val(R), Val(σ), VT)
-    stage_range!(py, px, pt, m, s, s - rem, s, 0, m, Val(R), Val(σ), Vec{2,T})
-  end
-end
-
-# Stage over a q range: the kernel behind the split of a single line, whose q extent is known
-# at generation time so the vectorised range and the scalar tail are emitted directly.
 # base + i·2^e as one scaled-index address computation (a GEP over an element of 2^e bytes), so
 # that the scaling stays in the addressing mode instead of becoming a value of its own.
 for (e, ty) in enumerate(("i8", "i16", "i32", "i64"))
@@ -89,6 +72,9 @@ stream_bases(base, r, R, first) =
 stream_strides(r, R, unit) =
   vcat(:($(Symbol(r, 1)) = $unit), [:($(Symbol(r, o)) = $o * $(Symbol(r, 1))) for o in (3, 5, 7) if o < R])
 
+# Stage over a q range: a stride that is not a multiple of the complex lanes of VT runs the
+# first s - s % C lanes in vectors of C complex values and the remainder one complex value at a
+# time.
 function stage_range_body(::Type{T}, R::Int, σ::Int, ::Type{VT}) where {T,VT}
   stmts, a, X = codelet(R, σ)
   C = ncomplex(VT); B = 2sizeof(T)
@@ -147,20 +133,6 @@ end
   stage_last_body(T, R, σ, VT)
 end
 
-# Last stage for strides that are not a multiple of the complex lanes of VT (odd batch sizes):
-# the first s - s % C lanes in vectors, the scalar remainder one complex value at a time.
-@inline function stage_last_split!(py::Ptr{T}, px::Ptr{T}, s::Int,
-                                   ::Val{R}, ::Val{σ}, ::Type{VT}) where {T,R,σ,VT}
-  C = ncomplex(VT)
-  rem = s % C
-  if rem == 0
-    stage_last!(py, px, s, Val(R), Val(σ), VT)
-  else
-    stage_last!(py, px, s, 0, s - rem, Val(R), Val(σ), VT)
-    stage_last!(py, px, s, s - rem, s, Val(R), Val(σ), Vec{2,T})
-  end
-end
-
 # First stage with s = 1 for vectors of C = 2 or 4 complex values (m a multiple of C):
 # butterflies p, …, p+C-1 share a register. Groups of C outputs X_k…X_(k+C-1) are transposed
 # (complex value l of each output into one vector) so that each store writes C consecutive