gelu: dispatch f16 through linalg and back it with a lookup table

[czoli1976]

GeluApproximate's f16 arm ran an inline scalar loop instead of calling ops().gelu_f16, so the f16 GELU dispatch slot was dead on every architecture — x86 has had x86_64_avx512_gelu_f16_16n registered but unreachable. This routes the canonical pow=3 f16 path through the dispatcher and backs it with a 2^16-entry lookup table.

Why a table, and why it is safe

Every f16 bit pattern is a valid index, so the whole activation becomes one load per element. The table is built from the same extracted gelu() scalar that HGelu8 evaluates, so it is bit-identical to the kernel it replaces by construction rather than by approximation — lut_matches_scalar_kernel_on_every_f16 checks that over all 65536 patterns.

It is built lazily behind a OnceLock: a model with no f16 GELU never allocates the 128 KiB.

A second test, registered_kernel_tracks_the_scalar_kernel_on_every_f16, holds whatever kernel is registered on the current arch to within 1 ULP of the scalar reference, so a per-arch override is checked on its own hardware rather than assumed.

Numbers

Kernel, cargo bench -p tract-linalg --bench gelu -- gelu_f16 (M-series):

buffer	scalar	lut
2 KiB	332 Melem/s	3.68 Gelem/s (11.1x)
128 KiB	329 Melem/s	3.56 Gelem/s (10.8x)
2 MiB	329 Melem/s	3.13 Gelem/s (9.5x)

It does not fall off a cliff when the 128 KiB table stops fitting alongside the data, so I also measured a transformer FFN block (matmul -> gelu -> matmul, f16) through a real plan, where the GEMMs compete for cache:

shape	scalar	lut
seq=128 d=768 ff=3072	3.12 ms	1.95 ms (1.60x)
seq=384 d=1024 ff=4096	17.89 ms	14.14 ms (1.25x)

First run is faster too (5.09 ms -> 2.47 ms), so the one-time table build repays inside the first inference. That FFN-only graph overstates GELU's share of a whole network — treat it as an upper bound.

Behaviour change, stated plainly

Two deltas, both from the rewire rather than the table:

• The old core scalar computed (2.0 / f32::consts::PI).sqrt() at runtime (0x3f4c4229); linalg uses the correctly-rounded literal (0x3f4c422a), 1 ULP apart. Across all 63488 finite f16 inputs this changes exactly one output, by one f16 ULP, and toward the true function.
• On x86 with AVX-512, x86_64_avx512_gelu_f16_16n now actually runs. It was already written and frame-tested, just never reached. I have no AVX-512 hardware, so the 1-ULP test above is what covers it — please look at that CI result specifically.

I have not run the large_models RBO nightly. Given #2318 was reverted for exactly this class of shift, that is the check I would want before this merges, and I could not run it here. Happy for it to sit until that is green.

Tests

tract-linalg 4410, tract-core 269, test-f16 2378, test-unit-core 816, plus nnef / hir / onnx / onnx-opl / transformers / extra / pulse — all green. cargo fmt --all clean; cargo clippy clean for the touched crates (the 3 --all-targets failures are pre-existing x86-only benches that do not build on aarch64).

Unrelated find

While benchmarking an alternative I measured arm64simd_gelu_f32_4n_fused returning 0.5 * x * 2^-23 instead of 0 for x below about -9 — its Pade tanh saturates one ULP short of -1, so 1 + tanh never cancels and the error grows linearly with |x| (-5.96e-5 at x=-1000, -3.90e-3 at x=-65504). That is pre-existing and affects f32, so it is not in this PR; I will open a separate issue with the clamp fix.

🍍