Performance and Numerical Reliability of Floating-Point Computation

Summary

Floating-point computation is fundamental to scientific computing, numerical simulation, and machine learning. Achieving high computational performance while maintaining numerical reliability remains a central challenge. Finite precision introduces rounding and approximation errors, while computational efficiency depends on numerical algorithms, implementation strategies, and hardware architectures. Understanding how these factors interact is essential for developing reliable and efficient numerical software.

Mathematical functions, such as exponentials, logarithms, and trigonometric functions, provide an important example of these challenges. Their implementations rely on numerical approximations and optimization techniques that introduce trade-offs between accuracy and performance. SIMD vectorization improves throughput by evaluating multiple inputs simultaneously, but the resulting numerical behavior and execution cost can depend on the function, the input values, and sometimes the other values in the same vector. Consequently, conventional benchmarks and documented accuracy bounds may not fully characterize an implementation’s behavior.

To investigate these effects, we systematically evaluate unary single-precision (FP32) mathematical functions in 256-bit AVX2 implementations from AOCL-LibM, GNU libmvec, SLEEF, and Intel SVML. For numerical accuracy, we exhaustively test every finite FP32 encoding against a 256-bit MPFR reference, supplemented by checks for infinities, NaNs, and mixed-lane inputs. We quantify errors in units in the last place (ULP) and examine whether implementations satisfy their documented accuracy guarantees. Selected tests also reveal accuracy coupling between SIMD lanes, where changing a neighboring input can alter the result for an otherwise unchanged value.

For performance, we evaluate inputs that exercise ordinary SIMD evaluation, expensive argument reduction, and scalar correction or fallback paths. Under identical benchmark conditions, every tested function-library pair achieves higher throughput through its vector interface on selected inexpensive inputs. However, input-dependent execution paths can reduce throughput by an order of magnitude or more, and some vectorized routines become slower than their scalar counterparts. We also identify violations of documented accuracy bounds in several tested implementations.

These findings demonstrate the importance of evaluating numerical accuracy and computational performance together, particularly across different input regions. Our results provide a reference for selecting and validating mathematical function implementations according to application requirements, while helping library developers identify numerical corner cases and performance bottlenecks. The measurements are specific to the tested library versions, FP32 AVX2 interfaces, and AMD EPYC system.

FP32 AVX2 accuracy

Libraries and interfaces

Tested 256-bit FP32 SIMD interfaces. * stands for the function name. The first six rows follow the paper’s library table; SVML CR is the later addition.

Backend Version AVX2 symbol family Error bound
AOCL 5.3.0 amd_vrs8_*f ≤4 ULP
libmvec glibc 2.42 _ZGVdN8v_*f ≤4 ULP (historical bound)
SLEEF U10 3.9.0 Sleef_*f8_u10avx2 ≤1 ULP (function-specific exceptions below)
SLEEF U35 3.9.0 Sleef_*f8_u35avx2 ≤3.5 ULP
SVML 2026.0.0-947 __svml_*f8 ≤4 ULP
SVML HA 2026.0.0-947 __svml_*f8_ha ≤1 ULP (conjecture)
SVML CR 2026.2.0-863 __svml_*f8_cr ≤0.5 ULP (conjecture)

The libmvec comparison uses the bound documented through glibc 2.41; glibc 2.42 removed that bound from its documentation. SLEEF bounds apply within each function’s documented effective region; U05, U15 and U35 exceptions are marked in the results. The SVML CR and HA bounds are conjuecture based on the compiler precision mode setting.

Maximum observed error

Observed maximum error in ULP within each implementation’s effective region; lower is better. Finite values are rounded upward.

Function AOCL libmvec SLEEF U10 SLEEF U35 SVML SVML HA SVML CR
acos 1.534 2.153 0.695 1.282 1.245 0.796 0.500
acosh — 2.807 0.513† — 2.807 0.902 —
acospi — — — — 2.457 0.792 —
asin 2.232 2.662 0.777 2.432 2.266 0.755 0.500
asinh — 2.087 1.001† — 2.087 0.716 —
asinpi — — — — 2.677 0.788 —
atan 2.504 2.466 0.939 2.637 2.466 0.829 0.500
atanh — 2.044 1.001 — 2.044 0.572 —
atanpi — — — — 3.269 1.000 —
cbrt — 2.852 0.501 1.312 2.774 0.522 0.500
cos 12.778 2.014 0.937 2.494 2.468 0.875 0.500
cosd — — — — 2.026 0.501 —
cosh 3.695 2.080 0.503† 1.447† 2.080 0.572 0.500
cospi — — 7.137e44 (U05) — 2.455 0.611 —
erf ∞* 1.929 1.000 — 1.929 0.780 —
erfc 6.628 1.615 1.448 (U15) — 1.615 0.934 —
exp 3.281 2.638 0.988 — 2.638 0.750 0.500
exp10 — 1.735 0.822 1.004 1.735 0.525 0.500
exp2 1.217 0.823 0.885 0.885 2.844 0.617 0.500
exp2m1 — — — — — — 0.500
expm1 — 1.679 1.000 — 2.620 0.560 0.500
invsqrt — — — — 2.386‡ 0.510‡ 0.500‡
lgamma — — 9.555e4† — — — —
log 1.999e6† 3.938 0.629 2.845 3.938 0.932 0.500
log10 4.635 2.785 0.712 — 2.785 0.771 0.500
log1p — 1.958 1.000 — 1.958 0.623 0.500
log2 4.295 2.082 0.674 3.073 2.082 0.681 0.500
log2p1 — — — — — — 0.500
sin 3.114 2.112 0.904 2.474 2.371 0.771 0.500
sind — — — — 1.959 0.501 —
sinh — 1.546 1.000† 3.403† 1.546 0.697 0.500
sinpi — — 8.389e6 (U05) — 2.454 0.611 —
sqrt 0.500 0.500 0.500 (U35) 0.500 0.500 0.500 —
tan 3.989 3.651 0.938 3.456 3.651 0.503 0.500
tand — — — — 3.913 0.880 —
tanh 0.500 1.027 1.001 2.860 3.819 0.514 0.500
tanpi — — — — 3.176 1.000 —
tgamma — — 2.463e4 — — — —

Reading the table

  • —: no result in the respective tested interface/campaign. exp2m1 and log2p1 are additions from the CR campaign.
  • †: a cutoff or special input increases the in-domain maximum beyond the effective-region value shown here.
  • ‡: invsqrt(-0) differs between SVML’s −∞ and the MPFR oracle’s +∞ convention. The CR scan records this one semantic discrepancy separately; it is not included in the finite effective-region maximum.
  • ∞*: AOCL erf can return infinity for a finite result with a NaN companion. Its maximum without NaN companions is 4.015 ULP; the broadcast maximum is 3.939 ULP. The infinity is a separate diagnostic, not the broadcast error.
  • U05/U15/U35: actual SLEEF routine suffix where it differs from the column heading.
  • libmvec sqrt uses a local AVX2 square-root intrinsic, because no native symbol is provided.

Code

The source code for our SIMD mathematical function evaluation is available on GitHub.

Publications

    Current people

    Contact