FMA dispatch and exp

The classical pair kernels return the same bits on every x86-64 instruction set, on AArch64, and with any C library. This page lists what that rests on.

Fused multiply-add

The pair-list force updates and the minimum-image fold use a fused multiply-add. A build for the x86-64 baseline picks the path once per process:

  • the hardware instruction, in functions compiled with target("fma"), when simd::detail::fma_hw_active() is true;

  • otherwise simd::exact_fma, a Dekker product and a round-to-odd sum (Boldo and Melquiond, IEEE Trans. Comput. 57(4), 2008).

fma_hw_active() reads CPUID leaf 1 (FMA, OSXSAVE and AVX) and XCR0 itself, so it needs no support from the C library. On glibc 2.33 and newer it also asks glibc, so GLIBC_TUNABLES=glibc.cpu.hwcaps-FMA= selects the software path. RGPOT_FMA=software selects it on any C library.

Both paths compute a*b+c with one rounding, so they return the same bits. ExactFmaEq checks this against the hardware instruction on 1.5 million cases: edge values, ties, near cancellation, the whole exponent range, and a million uniform bit patterns. FmaPathEq hashes LJ and Morse energies and forces on both paths, the software one through RGPOT_FMA=software.

exp

simd::exp_fixed replaces the C library’s exp in the Morse kernel. The C library’s exp has variants for each instruction set, and its last bit can change between them and between library versions. exp_fixed uses only IEEE additions, multiplications and bit operations, with no table and no fused multiply-add:

  1. Clamp x to [-746, 710], keeping overflow and underflow.

  2. Reduce x = k ln2 + r by the Cody-Waite split of fdlibm. k ln2_hi is exact and the subtraction is exact by Sterbenz’s lemma, so |r| < 0.3466=.

  3. Evaluate 1 + r + r^2 q(r), with q a degree 9 Horner polynomial.

  4. Scale by 2kas two powers of two, so only the last product rounds.

The scalar function and the vector loops run this body, and their bits match. exp_fixed_n and the Morse row kernel are compiled as target_clones for the baseline, AVX2 and AVX-512F. The loader picks one clone from libgcc’s own CPU check. Contraction is off in these files, and -fno-trapping-math -fno-math-errno let GCC vectorise them without changing a value.

Error bounds:

quantity

bound

source

polynomial against exp on the interval

3.8e-18 relative

Sollya supnorm

reduction, ln2 split

1.2e-26 absolute

exp_constants.py (sympy)

normal results, a priori, first order

1.24 ulp

exp_bound.py

measured against expl, 220inputs

0.944 ulp

ExpFixedEq on x86-64

The scripts are in scripts/exp_fixed/. fit.sollya reproduces the coefficients and the approximation bound. exp_constants.py checks the constants, the clamp limits and the exponent construction with exact rationals. exp_bound.py propagates a running rounding-error bound over the reduced argument. ExpFixedEq pins a digest of 220results, so a change in any bit fails the test on every platform.

The Morse kernel version moved to 2 with this function, so cached results computed with the C library’s exp are not reused.