Module fastmath 

Transcendentals without libm calls — the shared kernel every engine uses.

Why this module exists

exp, ln, tanh, sin, cos and powf have no SIMD instruction at any width. Every call is a scalar libm call, and — worse than its own cost — it is a hard barrier to vectorising the loop it sits in. Replacing them with polynomials is the highest-yield mechanical change available in this workspace (docs/plans/turbocharger.md).

sqrt, min, max, mul, add are the opposite: SSE2 baseline, already vectorised by the compiler. Rewriting those is how a campaign wastes a week — mp3’s xrpow hand-AVX2 measured 0.97x and was reverted.

Why ONE module and not three

Three implementations of this idea already existed, all wired to production, none aware of the others: ffai-diana’s exp_fast, ffai-mercury’s fast_exp, ffai-argus’s exp_poly. They had drifted in the one detail that decides whether the win happens at all (see below), and two of the three had the bug. Consolidating is the point, not a tidy-up.

★ The rounding step is the whole trick

exp(x) = 2^(x*log2 e) splits into an integer power (written straight into the f32 exponent field) and a fractional part (a degree-5 polynomial). That split needs a round-to-integer — and that is where two of the three implementations reintroduced the libm call they had just removed.

Rust’s f32::round is ties-away-from-zero, which no x86 instruction implements (vroundps is ties-to-even). So it lowers to a call, or to a long branchy sequence, sitting in the middle of the loop. floor is no better: it needs SSE4.1, above the portable x86-64 baseline.

Adding 1.5 * 2^23 forces any value below 2^22 to round into the mantissa’s last bit; subtracting it back leaves the value rounded to nearest-even. Pure float arithmetic — no call, no branch, and it needs no target_feature, so it vectorises on aarch64 (NEON is baseline) exactly as it does on x86.

Measured by ffai-diana over 16 M elements, best of 7, single thread:

	time	rate
memcpy (the roofline)	5.58 ms	24.04 GB/s
with f32::round	60.84 ms	2.21 GB/s
with round_ties_even	38.85 ms	3.45 GB/s
with the magic number	12.91 ms	10.40 GB/s

4.71x, and bit-identical to the round() version over the activation range. A transcendental within 2.3x of pure memory traffic is a transcendental that vectorised.

Accuracy and how it is gated

These are float approximations: the gate is a tolerance against libm plus the caller’s own end-to-end oracle, never bit-identity against std. Each function documents its measured worst case. The tests here check the tolerance, the landmarks, and the shape (monotonicity where it holds, saturation, exact values at 0) — the last of which catches an approximation that is accurate on average and wrong somewhere specific.

Functions

erf

erf(x), Abramowitz & Stegun 7.1.26 — 1.6e-6 measured in f32.

exp

e^x, accurate to 4.2e-6 relative over [-20, 20] (measured, not claimed — the degree-5 polynomial is the limit, and the two implementations this replaced both documented ~1e-7, which was optimistic).

exp2

2^x.

exp_sub_sum_inplace

row[i] = exp(row[i] - max), returning the sum — the softmax inner loop.

exp_sub_sum_scalar

The oracle. Every vector twin is gated against this.

gelu_erf

gelu in its exact form — 0.5x(1 + erf(x/sqrt 2)).

gelu_tanh

gelu_pytorch_tanh — 0.5x(1 + tanh(sqrt(2/pi)(x + 0.044715 x^3))).

gelu_tanh_inplace

xs[i] = gelu_tanh(xs[i]), vectorised.

ln

Natural log, accurate to ~1e-6 absolute over the positive range.

log10

log10(x), for log-mel and decibel work.

max_f32

Maximum of a slice — vectorised.

max_f32_scalar

The oracle.

round_ties_even_fast

Round to nearest even, without a libm call or an SSE4.1 instruction.

sigmoid

1 / (1 + e^-x).

silu

SiLU/swish — x * sigmoid(x).

tanh

tanh(x), via 1 - 2/(e^{2x} + 1).