How to Vectorize C++20 Transcendental Std Math Functions?

[ax3l]

Hi @mattkretz,

Our physics modeling in ImpactX uses a lot of trigonometric and generally transcendental functions like sin, cos, sinh, cosh, exp, pow, etc.

What is the status from your end for vectorizing them in vir-simd / ISO C++?
I read online and tested that auto-vectorizers are very good at this now, but I see no vectorization being picked up for them for my std::simd/vir-simd code written, and fallbacks to scalar ops.

Being so far only able to use C++20, I am exploring drafts (hacks?) like this (#53) right now (used here and here), but I was of course hoping vir-simd can solve this for physics-heavy code? :) Even on AVX2, this gives me a 2.37x speed-up in real-world use cases that we need, instead of mid-SIMD path unwrapping and sequentially solving the heavy math I came to use SIMD for.

I was also made aware that C++26 has vectorized these math functions now (?):

• https://eel.is/c++draft/simd#math
• https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2024/p1928r15.pdf

but I will need to get this into our C++20 code base and can maybe in 2-4 years switch to C++26 compilers due to HPC systems and large community stacks being slower (HPC systems in DOE, Helmholtz, Conda-Forge, etc.). What's the best way to pull in the vectorized transcendental math for now?

[ax3l]

I went ahead and prototyped this with Claude: #53. Short version is that the auto-vectorizer route works
but is not worth its side effects, and calling libmvec's entry points directly is.

Why the auto-vectorizer route was abandoned

Shaping each function as a short constant-trip lane loop does get GCC to collapse it into one
libmvec call, but only if two conditions hold, and both cost more than they look:

1. -fno-math-errno must be on the command line of the translation unit. A function-local
   #pragma GCC optimize("no-math-errno") does not work: the loop is vectorized in the caller's
   context after inlining, so the caller's flags decide.
2. The vector variants must be declared. glibc declares them only under -ffast-math, so the
   header has to declare them itself (__attribute__((__simd__("notinbranch"))), GCC only —
   clang ignores it and needs -fveclib=libmvec instead).

The killer is that both reach whole translation units. With the declarations in scope and the
flag set, GCC also auto-vectorizes ordinary scalar loops over sin/sinh into libmvec calls,
silently changing their accuracy. A full AMReX ctest sweep caught this as a failure in an
unrelated reduction test: two formulations of the same scalar sum, on a plain weighted moment
with no math function anywhere near it, drifted apart by one ULP because -fno-math-errno
changed what the compiler was willing to inline and contract.

Calling the vector entry points by name has none of this. No flags, no auto-vectorization, and
the effect is confined to the functions that opt in.

Why vir::vecmath and not vir::stdx

Declaring them in vir::stdx is tempting, because qualified lookup stops at the first namespace
with a match, so vir::stdx::sinh would resolve to the fast one with no ambiguity. I tried that
first. Three problems, all silent:

• It hides the whole overload set of that name. _GLIBCXX_SIMD_MATH_CALL2_ generates two
  overloads per two-argument function, the first with a second parameter deliberately excluded
  from deduction. Miss that and stdx::pow(x, 2.0) and stdx::pow(2.0, x) simply stop compiling.
• using namespace vir::stdx; becomes ambiguous for all 22 one-argument names, because the
  using-directive flattens both sets into one scope. Worse for pow/atan2: they resolve to the
  underlying overload and quietly take the slow path.
• Include-order ODR. An inline function calling vir::stdx::sinh gets a different body
  depending on whether its TU happened to include the header, under the same mangled name. The
  linker keeps one, so a TU that never mentions the header can change results with link order.

A namespace of its own has none of these, and vir::stdx is left exactly as it was. Note that an
unqualified sinh(x) still reaches the underlying implementation through ADL and no
using-declaration changes that, so callers have to name vir::vecmath::sinh. That is also why the
consuming side (AMReX) keeps its own amrex::Math:: layer as the single call site convention.

Portability

toolchain	result
GCC + glibc ≥ 2.35, x86-64	all 24 functions routed
GCC + glibc 2.22–2.34	sin cos exp log pow; the rest forward per lane
glibc < 2.22, non-glibc, non-x86-64, x32	inert, forwards per lane
vir's own simd (no <experimental/simd>)	inert
VIR_DISABLE_SIMD_VECMATH	inert

vir::vecmath declares all 24 names in every configuration; whether a given one is routed is the
header's business, not the caller's. The relevant glibc is the build-time one, so a toolchain
with an old sysroot (conda-forge ships 2.17) forwards even on a recent host — deliberately, since
unguarded it is an undefined reference to _ZGVdN4v_sinh at link time.

ISA selection was checked on SSE2, AVX-without-AVX2, AVX2 and AVX-512 targets: each picks the
right ISA letter, splits wider simds into whole chunks, and emits no -Wpsabi. The chunk callers
live in an ISA-named inline namespace so they do not collide across -mavx/-mavx2 TUs, and the
public overloads carry an ISA tag in the signature for the same reason.

One caveat worth knowing: _ZGVdN4v_* are ifuncs, so masking AVX2 at runtime
(GLIBC_TUNABLES=glibc.cpu.hwcaps=-AVX2) resolves them to SSE wrappers and makes sin slower
than libstdc++'s own. That bit us in CI before we spotted it.

Accuracy, errno and exception flags

libmvec is built for -ffast-math callers and glibc's own test suite exercises it with errno and
exception checking switched off, so both become unspecified. Measured against scalar libm over
2^18 points per function: ≤ 3 ULP everywhere (glibc documents 4), sqrt exact, and zero
sign-of-zero, ±Inf or NaN mismatches. Denormals are not flushed.

What does differ: sin(inf) sets errno where the scalar routine does not, exact results such as
cbrt(-8) pick up FE_INEXACT, atanh(1) reports FE_INVALID for what is a pole error, and
because flags are per call rather than per lane, one lane hitting a pole leaves them set for the
whole vector.

This lines up with where the standard is heading. P1928R15 §6.1: "The intent is to avoid errno
altogether, while still supporting floating-point exceptions (possibly depending on compiler
flags)", noted as needing more work and not yet in the wording. So dropping errno is the
direction of travel; the exception flags are where a vector math library falls short of it. No
accuracy bound is specified at all, so the ULP figures are glibc's documentation, not a guarantee.
That same section also says the basic_simd overload wording "leaves room for ambiguity", which
matches what the pow/atan2 overload set did to me above.

Test coverage

Two files, registered in ext_tests rather than simd_tests — the latter only runs when vir uses
its own simd, which is precisely the configuration where this header is inert.

simd_vecmath.cc checks all 24 functions against the scalar routines for every element type and
ABI the harness instantiates, with the tolerance chosen from whether the call is actually routed:
4 ULP where libmvec answers, 1 where the underlying implementation does. The widths that decide
chunking (1…32) are swept separately, lane by lane, because the harness's ULP helper does not
accept fixed_size ABIs on top of libstdc++'s simd. Also covered: all four spellings of the
two-argument functions, the names the header deliberately leaves alone (hypot, sqrt, abs),
taking a function's address, and the chunk-width selection as compile-time assertions.

simd_vecmath_disabled.cc asserts VIR_DISABLE_SIMD_VECMATH really disables it, by identity
against the underlying implementation rather than against scalar libm — comparing against libm
there would measure libstdc++'s argument reduction, which returns inf for cos(finite_max).

Mutation-checked: truncating the chunk loop, restoring an unqualified self-call in the reversed
two-argument overload, removing the fallback overload, and narrowing the chunk width are each
caught. The first two were real bugs the tests found while being written — the reversed overload
in particular compiled and returned correct answers while silently taking the slow path.

Not coverable in the suite and verified by hand: the multi-TU ODR case (needs two objects built
with different -m flags and a link) and AVX-512 execution (no hardware here, checked by reading
the generated assembly).