Add an estimating planner to the NEON, SSE and wasm_simd planners

[HEnquist]

Draft, and stacked on #179, which is also still a draft. A PR here cannot be based on a branch in my
fork, so the diff below contains #179's commits as well. The estimating planner starts at ba8d179,
"Bring over the planner tuning harness", and is about 4700 of the 7000 lines.

On this description: I asked Claude for a detailed description of how the planner works, and it
went all in. So this is a lot more text than I would have written myself. The numbers are real and
reproducible with the commands under Results. If the length puts you off, read the TL;DR below, and
then "Where it loses" and "Known gaps" if you want the caveats. If there is a number you want that
is not here, just ask and I will measure it.

TL;DR

1. Enumerate the recipes that could compute a length, price each one, keep the cheapest. NEON, SSE
   and wasm_simd only; the scalar and AVX planners are untouched.
2. NEON f32, every length from 1 to 1000: faster at 471 of them by more than 5%, and at 208 by
   more than 20%, for a geomean of 0.891. 33 lengths lose by more than 5%. The same range in f64
   wins at 343 and loses at 2.
3. Faster on average in all six machine and element type combinations measured, at a geomean of
   0.88 to 0.98 of the fixed planner's runtime. Sampling 334 lengths up to a million instead, it
   wins by more than 5% at 64 to 161 of them and loses by more than 5% at 3 to 20.
4. Best case 5.46x faster (135721 on a Pi 5, f64). Worst case 3.51x slower (774209 on a Pi 5, f32),
   which is the one loss mode that needs the cache size to fix properly.
5. Planning costs 113x what the fixed planner's does, which is 4.3% of plan-plus-build.

What it replaces is a fixed decision tree. That fixed planner stays in the tree behind the
non-default tuning feature so the two can be measured against each other, which is what every
number here is. The losses fall into three modes, all described below, and none is fixed in this PR.

Results

ratio throughout is the estimating planner's runtime over the fixed planner's, so below 1 is
faster. A win or loss means more than 5% in that direction. Machines: Apple M1 (12 MB L2),
Raspberry Pi 5 (2 MB L3), i3-8100T (6 MB L3). Repeat runs hold the geomeans to about 0.001 but move
the win and loss counts by one or two, so read any count near the 5% threshold as approximate.

Survey: 334 lengths drawn from 1 to 1,000,000

survey --count 350 --seed 1 1..1000000, the same lengths on every machine.

backend / machine	geomean	p10	p50	p90	worst	wins	losses	same recipe
NEON f64, M1	0.917	0.76	0.98	1.000	1.14	147	5	120
NEON f32, M1	0.933	0.78	1.00	1.000	1.40	118	15	141
NEON f64, Pi 5	0.881	0.70	0.96	1.000	1.26	161	3	120
NEON f32, Pi 5	0.938	0.79	1.00	1.000	3.51	123	14	141
SSE f64, i3-8100T	0.983	0.91	1.00	1.016	1.41	64	20	175
SSE f32, i3-8100T	0.943	0.81	1.00	1.000	1.38	110	7	158

SSE f64 is the weakest cell and the only one whose p90 is above 1. It agrees with the fixed planner
at 175 of 334 lengths, wins at 64 and loses at 20.

Dense sweep: every length from 1 to 1000, on the M1

The survey samples, so it can miss a whole neighbourhood. This trades range for completeness: a much
narrower band, but every length in it. Note that no length here comes near cache_elems, so loss
mode 1 below cannot appear at all, and these rows are not comparable in difficulty with the survey.

	lengths	same recipe	geomean	wins >5%	wins >20%	losses >5%	worst
NEON f64	1000	533	0.934	343	125	2	1.09
NEON f32	1000	308	0.891	471	208	33	1.21

Where the wins come from

Three recipe changes account for most of them.

1. Rader's instead of Bluestein's at a prime. Now that the Rader's permutation is a table
   (Precompute the Rader's permutation #178), Rader's is the better answer far more often than the fixed planner assumes. Length 523
   goes from bs(523,r4(3,b24)) to rad(rn(6.3,b29)), 2.03x faster on the M1 in f64. It compounds:
   933 and 622 win 1.72x and 1.71x purely because their inner 311 changed.
2. A better RadixN base. rn(7,mrs(b3,b9)) becomes rn(7.3,b9) at 189, 1.80x. In f32 this
   class is large: 448 goes from rn(2.4.4,gts(b2,b7)) to rn(7.2,b32), 2.07x.
3. Radix4 instead of a butterfly pair. 48 goes from mrs(b6,b8) to r4(1,b12), 1.66x in f64
   and 2.05x in f32, and 96 from mrs(b8,b12) to r4(1,b24), 1.97x in f32.

The largest single win is 135721 on the Pi 5 in f64, 5.46x, where bs(135721,r4(7,b24)) becomes
rad(rn(6.6.5.2,gts(b13,b29))). The same length wins 2.32x on the M1 and 1.78x on the i3, which is
the shape of the whole Bluestein's-to-Rader's class: biggest where cache is smallest.

Where it loses

Broken out by the largest prime factor of the length, which is what separates the loss classes:

class	M1 f64	M1 f32	Pi 5 f64	Pi 5 f32	i3 f64	i3 f32
every factor has a butterfly	0.950 / 1	0.937 / 2	0.950 / 2	0.941 / 2	0.977 / 3	0.956 / 0
largest prime 33 to 2000	0.925 / 0	0.912 / 7	0.906 / 0	0.918 / 2	0.983 / 7	0.904 / 2
largest prime 2k to 20k	0.873 / 2	0.954 / 1	0.852 / 1	0.988 / 3	1.002 / 8	0.971 / 1
largest prime above 20k	0.835 / 2	0.961 / 5	0.645 / 0	0.936 / 7	0.978 / 2	0.998 / 4

Each cell is geomean / losses beyond 5%. There are three distinct loss modes.

1. One large Bluestein's whose inner FFT no longer fits in cache. This is the worst class and
every worst case in the table above is in it. The model takes a single Bluestein's over the whole
length where the fixed planner splits first and keeps its inner FFTs cache resident:

len	machine	ratio	fixed	model
774209	Pi 5 f32	3.51	mr(rad(rn(5.3,gts(b2,b11))),bs(2339,r4(4,b24)))	bs(774209,r4(8,b24))
232371	Pi 5 f32	2.58	mr(bs(25819,r4(6,b16)),b9)	bs(232371,r4(7,b32))
232371	i3 f64	1.41	rn(3.3,bs(25819,r4(6,b16)))	bs(232371,r4(7,b32))
969365	i3 f32	1.38	mr(b5,bs(193873,r4(7,b24)))	bs(969365,r4(8,b32))
616197	M1 f32	1.17	mr(bs(205399,r4(7,b32)),b3)	bs(616197,r4(8,b24))

At 774209 the model's inner FFT is 1572864 points, 12.6 MB in f32. Timing the pieces on the Pi puts
out-of-cache work at 0.79 to 0.98 ns per unit of cost against 0.19 to 0.21 for cache-resident work,
so it is still under-priced about four times over after dram_pass. The same pick is correct on
the M1, where 12.6 MB fits in a 12 MB L2. The multiplier is not the machine-specific part, the
threshold is, and cache_elems is a compile-time constant here. See "Known gaps".

2. The general-versus-small transpose crossover, in f32. 29 of the 33 losses in the f32 dense
sweep are one substitution, the model taking the general MixedRadix or GoodThomas where the
fixed planner takes the Small variant, at lengths 195 to 992:

195  1.18x   gts(b13,b15) -> mr(b13,b15)
209  1.18x   gts(b11,b19) -> mr(b11,b19)
228  1.21x   gts(b12,b19) -> mr(b12,b19)
234  1.20x   rn(3.3,gts(b2,b13)) -> mr(b13,rn(3,b6))

general_row was fitted against twelve measured general-over-small ratios in f64 and put the
crossover near length 200. In f64 this mode produces no losses at all over the same 1000 lengths;
in f32 the crossover clearly belongs higher. It is one weight, and making it per element type is
the obvious fix, but it wants its own measurement pass and is not in this PR.

3. Very short lengths where RadixN replaces a butterfly pair. The only two f64 sweep losses:
21 at 1.09x (gts(b3,b7) to rn(3,b7)) and 20 at 1.06x (gts(b4,b5) to rn(2,b10)). radix_call
prices the generic driver's fixed overhead but evidently not quite high enough at the bottom end.
Both stay under 10%, and they are the entire f64 loss list for lengths 1 to 1000.

Plan time

Enumerate-and-price costs far more than the fixed planner's plan. That is the wrong denominator: a
caller pays plan plus build, and building dominates. On the M1 in f64:

    len   plan fixed  plan estimate        build   extra vs plan+build
     32          351            252          664         -9.8%
     64          212            466         1394         15.8%
     96          302           6805         1527        355.5%
    256          173            424         2734          8.6%
   1000          499          15972         8792        166.5%
   4096          166            390        26681          0.8%
  10007          362         169077       477133         35.3%
  65536          156            372       376725          0.1%
 100000          347          49589       570450          8.6%
1000000          260          75949      5854092          1.3%

  totals: plan fixed 2828 ns, plan estimate 319295 ns (112.9x), build 7320191 ns
  the estimating planner adds 4.32% to plan-plus-build

Planning alone is 113x dearer, and that is 4.3% of plan-plus-build over this set. The worst relative
case is a short length with many divisors (96, 1000), where the absolute cost is single-digit
microseconds. Butterfly lengths and powers of two short-circuit enumeration entirely, which is why
32, 64, 256, 4096 and 65536 stay flat. These are cold-start figures with a fresh planner per
length; a planner reused across lengths shares inner lengths through its caches.

Reproducing all of the above

cd tools/planner_tuning && cargo build --release
./target/release/planner_tuning survey --count 350 --seed 1 1..1000000
./target/release/planner_tuning survey --count 350 --seed 1 --f32 1..1000000
./target/release/planner_tuning sweep 1..1000
./target/release/planner_tuning sweep --f32 1..1000
./target/release/planner_tuning plantime 32 64 96 256 1000 4096 10007 65536 100000 1000000

The tool crate selects the SIMD feature for the host, neon on aarch64 and sse on x86-64. Do not
build the library with default features for this, or the planner would pick AVX and measure the
wrong backend.

How it works

Three steps, and only the middle one is new. All three are in design_fft_for_len in each planner.

1. Enumerate. List the recipes that could compute this length as Shape values: a top-level
   algorithm plus the lengths of its inner FFTs, but not the inner recipes.
2. Price. Give each candidate a number, with CostModel::cost.
3. Keep the minimum. Ties go to the first candidate, which is always the fixed planner's pick.

The model lives at src/simd/simd_estimate.rs and is shared by all three SIMD planners. The scalar
and AVX planners do not use it and are untouched.

What a cost is

One unit is one issued arithmetic instruction. A cost is not nanoseconds and nothing converts it
to time, because only the ranking within one length is ever used. That is why one weight set travels
across machines of different clock speeds.

cost = counted arithmetic instructions        (read off the source)
     + a memory term                          (accesses x pattern x whether it fits in cache)

A pure operation count, the FFTW_ESTIMATE analogue, scores worse than the fixed planner. The
memory term is what makes the model work; it is not a refinement.

The memory term

Each pass is charged per element it touches, times a multiplier for how it walks memory:

pattern	weight	what it is
sequential	1.0	a contiguous run, one cache line feeding many elements
strided	strided	a fixed stride, as in a transpose or a cross-FFT layer
permuted	permuted	a computed address per element: digit reversal, CRT, Rader's

Sequential and strided accesses are divided by the number of complex numbers in a vector; permuted
ones are not, because a gather computes an address per element and cannot fill a vector. Above
cache_elems complex numbers an access costs dram_pass extra, and a transpose's accesses cost
dram, because a cross layer keeps its locality inside the chunk it is already working on while a
transpose walks the whole rectangle.

What each node costs

node	arithmetic	memory
Butterfly(len)	counted table lookup	2*len sequential
Radix4 { k, base_len }	reps x base, plus per layer len/4 butterfly4 and 3 twiddles	2*len permuted digit reversal, 2*len strided per layer, plus radix_call
RadixN { factors, base_len }	reps x base, plus per layer len/r x (butterfly r + r-1 twiddles)	same as Radix4, plus radixn_extra per element per layer
MixedRadix	h x left + w x right, plus len twiddle multiplies	three transposes plus one sequential pass
GoodThomas	h x left + w x right, no twiddles	two permuted CRT passes plus one transpose
Raders { len }	2 x inner, len twiddles, len * rader_index	two permuted passes plus one sequential
Bluesteins { len, inner_len }	2 x inner, inner_len pointwise multiplies	sequential over the inner length and twice over the outer

Three terms exist only to encode a decision the model would otherwise be unable to make:

1. small versus general. The Small variants call transpose_small, a naive strided double
   loop, so their transposes are priced Permuted. The general variants hand the job to the
   transpose crate, which tiles the rectangle and pays general_row per row for it. That makes
   the general form cheaper per element and dearer per row, which is a crossover. Without the
   per-row term the model has only per-element costs and would take the general form at every
   length.
2. Which ordering of a split. small_row * max(width - height, 0) is the only thing separating
   mrs(A,B) from mrs(B,A). transpose_small's outer loop runs width times, so both Small
   variants change by exactly width - height when the pair is reversed. Charging the difference
   keeps a square pair free and leaves the better ordering at the cost it had before the term.
3. A generic driver against a hand-written one. radix_call is what a RadixN or Radix4
   execution costs regardless of length, and radixn_extra is what its cross layers cost per
   element over Radix4 doing the same work. Without the first, short lengths take a RadixN that
   measures 1.32x slower than a table-driven GoodThomasAlgorithmSmall.

Assumptions

These are the places where the planner asserts something rather than measuring it. Each was checked
against measurement, and the check is the thing to redo if a kernel changes.

1. A length with its own butterfly, or a power of two, has nothing to decide. Enumeration
   answers immediately at both. Over lengths 8 to 128 on NEON the bare butterfly is fastest at every
   one, and at every power of two from 64 up, across four datasets, the fixed planner's Radix4 is
   exactly the fastest measured candidate. This shortcut is also what keeps plan time flat at
   exactly the lengths where plan time is the largest fraction of plan-plus-build.
2. Only the smaller-width ordering of a two-way split is offered. The reverse roughly doubles
   the candidate count at a highly composite length for almost no information: the smaller-width
   ordering is the better one in 90 to 97% of measured pairs for the Small variants, and the
   general variants are usually indistinguishable.
3. Bluestein's is offered only where some prime factor has no butterfly of its own. If every
   factor has one, the whole length decomposes into butterflies, and across four sweeps of 1 to 1000
   not one such length was won by Bluestein's. It is offered at composite lengths: 671 = 11 x 61
   measured 1.46x faster as Bluestein's than as a split around a Rader's for 61.
4. A node's cost depends only on its own subtree. A node is charged at the size of the buffer it
   walks itself, never at the size of the transform it is nested inside. This is what makes the
   search affordable: each planner memoises the best cost per length beside its recipe cache and
   recurses over the divisors of a length rather than over whole trees. If a future term ever
   makes a node's cost depend on its parent, that memoisation becomes wrong, and the failure would
   be a subtly bad inner recipe rather than anything that trips a test.
5. Ties go to the fixed planner. Its pick is always the first candidate, so a change that does
   not move a cost cannot move a plan.
6. At most 48 candidates per length. Above that, only splits are dropped, most lopsided first,
   on the grounds that a split with a tiny side is mostly its large side plus a transpose.
   Everything structural is kept.
7. One weight set per backend and element type, not per machine. The weights are compile-time
   constants. Where a weight's optimum differs between machines it was swept on all of them and set
   to the value whose worst machine looks best.
8. wasm_simd borrows NEON's instruction counts. A quick test on an M1 showed wasm tracking NEON
   closely, but that is the same machine, and nothing has been counted or fitted for wasm itself.
   This is the weakest assumption in the PR and it is marked as a placeholder in the source.

What is tuned from measurements

The distinction the whole approach rests on. Most of the model is read off the source and must be
maintained when the source changes. A small set of weights is fitted against measurements and must
be rechecked when the machines change.

Read off the source, never fitted

quantity	where it is read
butterfly instruction counts	src/neon/*.rs, src/sse/*.rs, by hand; derived in OP-COUNTS.md
generated prime butterflies	the generator's loop structure, as a closed form that cannot drift
mul_complex, column_butterfly4	the vector trait impls
access counts and patterns per pass	the algorithm source, loads plus stores per pass

Fitted, and where each came from

weight	value	how it was chosen
strided	1.5, or 2.5 on SSE f32	grid against measured dumps
permuted	2.5	grid against measured dumps
general_row	30	the twelve measured general-over-small ratios
small_row	10	one outer iteration of transpose_small, 1.48 ns on an i3, 0.7 to 1.0 on an M1. A tie-break with a derivation; anything from 2 to 24 scores the same
rader_index	2 on NEON, 20 on SSE	one load from the permutation table #178 added, plus assembling the element, which costs far more on SSE. It was 30 and 45 when that index came from a loop-carried modular multiply
radixn_extra	0 on NEON, 6 and 1 on SSE f64 and f32	expected near zero where 2R rows fit the register file: 32 vector registers on aarch64, 16 on x86-64
radix_call	100	sized at length 14 on an M1, about 8 ns or 25 cycles: the call, the scratch split, the layer setup and the virtual call into the base FFT
cache_elems	256 KiB worth	the smallest last-level cache worth planning for, not any one machine's. No pick below length 16385 changes at this threshold, so weights fitted on short lengths stay valid
dram, dram_pass	6 and 2	fitted on a Pi 5, checked on an M1 and on an i3-8100T whose 6 MB L3 sits between the two

dram and dram_pass are the first weights whose optimum depends on the machine rather than the
backend, because they price the memory system rather than the instruction set. They have to move
together: charging ordinary passes without charging transposes just as hard makes a MixedRadix
wrapped around a smaller radix recipe look good, and those recipes measure worse on both machines.

What to redo when the kernels change

In this order, from free to expensive. The first four need no machine.

1. verify at a spread of lengths. It checks every enumerated candidate against an f64 reference
   DFT, which catches an illegal spec such as a Bluestein's inner shorter than 2n - 1. Never trust
   a timing taken before verify is clean.
2. explain on a recipe whose cost you can predict.
3. picks before and after, over tens of thousands of lengths in seconds. It says how far the
   change reaches before anything has been measured.
4. score against an existing dump, to see whether the ranking moved.
5. sweep 1..1000 and a survey up to a million if it did. Both cost-model defects found during
   development were invisible to the 33- and 44-length tuning sets and showed up only in a sweep.

A dump goes stale when the algorithm it measured changes. Every dump taken before #178 holds the
old, slow Rader's, so it cannot judge any Rader's decision.

Known gaps

1. A large prime factor, in f32, on a machine with a small cache, described under loss mode 1
   above. The fix, when it is wanted, is to read the last-level cache size at construction and set
   cache_elems from it: /sys/devices/system/cpu/cpu0/cache/ on Linux, and
   hw.perflevel0.l2cachesize on macOS, whose flat hw.*cachesize keys report the E-core sizes.
   At a 12 MiB threshold, dram_pass 3 costs the M1 one loss instead of seven while still fixing
   the Pi, so giving each machine its own cache size makes the conflict disappear.
2. The general-versus-small crossover is not per element type, which is loss mode 2 and 29 of
   the 33 f32 sweep losses.
3. Width and height are tied. mr(A,B) and mr(B,A) get the same cost apart from the
   small_row tie-break, yet many reversed pairs measure more than 2% apart. This is the clearest
   unexploited improvement and it is derivable from the code.
4. MixedRadix versus GoodThomas is a fixed preference, not a decision. With identical inners the
   cost difference is a constant multiple of len, so its sign cannot vary with length. A
   stride-aware rewrite aimed at this regressed both backends and was reverted.
5. x86 has one machine. The Pi 5 separated machine from backend on ARM; nothing has done that
   for SSE, so an SSE weight and an i3-8100T weight are still the same column. rader_index shows
   why that matters: 2 on NEON against 20 on SSE is a gap explained by how each backend gathers a
   complex number, but only a second x86 machine can confirm that reading.
6. wasm_simd has no counts or measurements of its own, per assumption 8.

What is in the diff

src/simd/simd_estimate.rs	the shapes, the enumeration and the cost model, shared by the three SIMD planners
src/neon, src/sse, src/wasm_simd planners	the enumerate-price-pick loop, the cost cache, and set_estimating for comparison
src/tuning/	exhaustive enumeration and planner adapters, behind the non-default tuning feature, #[doc(hidden)]
tools/planner_tuning/	the measurement harness, not shipped, plus COST-MODEL.md, OP-COUNTS.md and README.md

The harness drives the planners through the tuning feature, so recipes are built by the planners'
own code and priced by the library's own cost model: what gets timed is exactly what a planner would
construct. dump times every candidate at a length and writes a TSV; score, costs and explain
then replay that file offline, so iterating on the cost model needs no machine after the first run.
The exhaustive enumeration in src/tuning/ is deliberately separate from the planner's, because
scoring a pick needs alternatives no planner would propose.

Open question

The fixed planner is still in the tree, reachable only from the tuning feature. If this lands it
should probably come out, but keeping it is what makes every number above reproducible, so I have
left that decision for review.