Optimization: SIMD polynomial exp (float, 8-wide) + FTZ/DAZ

Author: Linus
Predecessor: variant_float_inlining_SoA_AVX2.cpp
File(s): variants/variant_float_inlining_SoA_AVX2_polyexp.cpp, main-float.cpp (FTZ/DAZ enable)
Status: TBD, positive locally (approx twice as fast)

What this changes

The 8 scalar expf calls in the j-loop become a single SIMD polynomial approximation exp256_ps that operates on __m256. Range-reduce x = n*ln2 + r with n = round(x/ln2), evaluate a degree-5 Horner polynomial for exp(r) on r in [-ln2/2, ln2/2], then scale by 2^n via the Schraudolph bit-trick on the IEEE 754 exponent slot (shift the biased exponent by 23). Degree-5 is enough for float, the polynomial is good to ~1e-5 absolute on the reduced range, well under the 1cm verify tolerance.

main-float.cpp also gets two MXCSR writes at the top of main():

_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);

This was the part that actually mattered, see Analysis.

Code (relevant snippet)

See variants/variant_float_inlining_SoA_AVX2_polyexp.cpp for the full j-loop and main-float.cpp for the FTZ/DAZ setup.

Benchmark

Conditions: min of N runs (TBD), N_ITER=1000, -O3 -march=native -flto. Diff is vs variant_float_inlining_SoA_AVX2.

x86

Machine: TBD

Frequency: TBD

N	W(n) [Flops]	T(n) [Cycles]	Performance [Flops/Cycle]	Diff vs Predecessor
50
100
200
500

ARM (Apple)

Not applicable, AVX2-only. A NEON port would need its own vectorised exp.

Analysis

roughly twice as fast

Takeaway

Polyexp helps as expected. The real fix is FTZ/DAZ in main-float.cpp, without it everything downstream in the float chain is silently subnormal-stalled

References

Schraudolph, N. N. (1999). "A Fast, Compact Approximation of the Exponential Function." Neural Computation 11(4), 853-862.
IEEE 754 double-precision format. Wikipedia. https://en.wikipedia.org/wiki/Double-precision_floating-point_format
Horner's method. Wikipedia. https://en.wikipedia.org/wiki/Horner%27s_method

2.5 KiB Raw Blame History