[WIP] Add simd_vecmath.h: transcendentals via a vector math library - #53
Open
ax3l wants to merge 12 commits into
Open
[WIP] Add simd_vecmath.h: transcendentals via a vector math library#53ax3l wants to merge 12 commits into
simd_vecmath.h: transcendentals via a vector math library#53ax3l wants to merge 12 commits into
Conversation
SIMD hardware has instructions for sqrt and abs but not for the transcendental functions, so simd math falls back to calling the scalar libm routine once per lane. That makes a vectorized kernel calling sin, exp or sinh slower than the scalar one it replaced: measured at width 4 on an i9-12900H, stdx::sinh took 8.8 ms per 2^20 evaluations where a plain scalar loop took 8.0 ms. glibc ships libmvec, which computes a whole register per call. This header calls its entry points directly, by their x86-64 vector function ABI names, so no compiler flags and no auto-vectorization are involved. The same measurement drops to 1.2 ms, 7.2x faster than before and 6.4x faster than the scalar loop. The overloads go into vir::stdx, where qualified lookup finds them ahead of the underlying implementation without ambiguity, so vir::stdx::sinh(x) simply becomes fast. Because a declaration in vir::stdx hides the underlying ones completely, each function comes as a pair: the vector math path for the element types and widths libmvec covers, and a fallback that forwards everything else, which is what keeps long double, odd widths and the scalar ABI working. Covers the 22 one-argument functions plus pow and atan2. Not hypot: the underlying implementation already evaluates it with SIMD instructions, including fixups libmvec's version does not do, and its overload set is larger than the pair generated here. Verified on AVX2, AVX-without-AVX2, SSE2 and AVX-512 targets, each selecting the right ISA letter and splitting wider simds into whole chunks, with no spills beyond what a multi-chunk call needs. Accuracy is within 3 ULP of scalar libm, inside libmvec's documented 4 ULP. Compiles clean under C++17 and C++20, GCC and clang, -Wall -Wextra, and stays inert on toolchains without a suitable glibc. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A portability review of the previous commit turned up four defects, three of them consequences of the same thing: declaring a name in vir::stdx hides the underlying overloads completely, so the replacements have to match what they displace exactly. pow and atan2 lost their converting forms. libstdc++ generates two overloads per two-argument function, the first with a second parameter deliberately excluded from deduction, which is what lets pow(x, 2.0) broadcast the scalar, plus a reversed form for pow(2.0, x). Deducing both parameters compiled fine and silently took all four spellings away: stdx::pow(x, 2.0), stdx::pow(2.0, x), stdx::atan2(x, 1.0) and stdx::atan2(1.0, x) stopped compiling for anyone who included this header. Follow libstdc++'s shape instead. The note on floating point semantics was simply wrong. It claimed libmvec neither sets errno nor raises exceptions; on glibc 2.39 it does both, inconsistently. 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 errno set for the whole vector. Say what actually happens, and say it next to the accuracy caveat. Add VIR_DISABLE_SIMD_VECMATH. Without an opt-out, a user who needs errno, the exception flags, sub-ULP accuracy, or who lands on a glibc built --disable-mathvec, has no recourse once the header arrives transitively. Give the emitted symbol an ISA discriminator. These functions are always_inline, so no out-of-line copy exists at -O1 and above, but at -O0 or when the address is taken a weak symbol appears, and two translation units built for different instruction sets agreed on its name while disagreeing on its body. The linker kept one, leaving the other TU calling a libmvec entry point its target may not be able to execute. libstdc++ guards its own math functions the same way with __odr_helper. Also test __GLIBC__ only after <features.h> has been included rather than relying on simd.h having pulled in <cstdlib> first, and note that the AVX entry points are wrappers around the SSE2 routine rather than true 256-bit implementations. Re-verified: the four overload spellings compile again, symbols now differ per ISA, the opt-out is inert, and the SSE2, AVX, AVX2 and AVX-512 targets still select the right entry points with no warnings. Accuracy and performance are unchanged (kernel benchmark 3.6x over scalar, no compiler flags). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Two test files, registered as ext_tests so they run against libstdc++'s simd. The general simd_tests pool only runs when vir falls back to its own implementation, which is precisely the configuration where this header does nothing, so putting them there would have tested an empty header. 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 which path the simd takes: 4 ULP where libmvec answers, 1 where the underlying implementation does. Since the harness only instantiates the scalar and native ABIs unless the expensive tests are enabled, and its ULP helper does not work with fixed_size ABIs on top of libstdc++'s simd, the widths that decide chunking are swept separately, from 1 to 32, with a lane-wise comparison that needs nothing from the harness. Also covered: all four spellings of the two-argument functions, the names the header deliberately leaves alone, taking a function's address, and the chunk width selection. simd_vecmath_disabled.cc asserts VIR_DISABLE_SIMD_VECMATH really disables it, by identity against the underlying implementation rather than against the scalar routines. The distinction matters: libstdc++'s vector cos returns inf for finite_max, so a comparison against libm would have been measuring the underlying implementation rather than the opt-out. Writing the tests turned up two defects. The reversed two-argument overload called itself unqualified, so argument-dependent lookup added the underlying overload, which deduces both parameters and therefore wins partial ordering: pow(2.0, x) compiled, gave correct answers, and quietly took the slow path. And chunk width selection had no coverage at all, because narrowing a chunk stays correct and only costs calls; it is now four compile-time assertions, the load-bearing one being that a width the native register divides has to use the native chunk. Mutation-checked: truncating the chunk loop, restoring the unqualified call, removing the fallback overload and narrowing the chunk width are each caught. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A second review round found that declaring these in vir::stdx costs more than
it buys. Qualified lookup does find them there without ambiguity, as intended,
but that is not the only way the names get looked up:
- `using namespace vir::stdx;` made all 22 one-argument names ambiguous,
because the using-directive flattens both these overloads and the ones the
underlying implementation contributes into a single scope. pow and atan2
were worse than ambiguous: they resolved to the underlying implementation
and quietly took the slow path.
- An inline function calling vir::stdx::sin got a different body depending on
whether its translation unit happened to include this header, with the same
mangled name. The linker keeps one, so a translation unit that never
mentions the header could silently change results and errno behaviour with
link order.
- And because a declaration in vir::stdx hides the whole overload set of that
name, every form the underlying implementation offers had to be reproduced
exactly. Getting that wrong was silent, as the earlier pow(x, 2.0)
regression showed.
Move them to vir::vecmath, called qualified. None of the above applies there,
and vir::stdx is left exactly as it was. An unqualified sinh(x) still reaches
the underlying implementation through argument-dependent lookup, and no
using-declaration can change that, so callers name what they want.
vir::vecmath now offers the same names whatever the glibc: functions without a
vector variant on this system get a forwarding overload rather than not being
declared. Whether a function is routed is this header's business, not the
caller's.
Also from the review: the chunk callers live in an inline namespace named after
the instruction set, so they no longer collide across -mavx/-mavx2 the way the
public overloads no longer do; x32 is excluded, where __x86_64__ is defined but
the ABI is not this one.
The tests follow the namespace, and now stand down where the header is inert
rather than testing the underlying implementation in its place. That was
failing on every platform without libmvec: libstdc++'s vector cos returns inf
for finite_max, so comparing it against scalar libm measured argument reduction
rather than anything to do with this header. Both files also no longer name
exp, exp2, expm1, cbrt or fabs unconditionally, which vir's own simd does not
provide, so they compile where <experimental/simd> is absent. pow's spelling
checks use a non-integral exponent, since with 2 the vector and scalar results
agree bit for bit and the comparison proved nothing.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
P1928R15 section 6.1 says the intent is to avoid errno altogether while still supporting floating-point exceptions, and that this needs more work and is not in the wording yet. Worth citing next to the caveats here: dropping errno is where the standard is heading, the exception flags are where a vector math library falls short of it, and no accuracy bound is specified at all, so the ULP figures are glibc's own documentation rather than a guarantee. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Not up to standards ⛔
|
Adds a user-facing section: what it does, which functions it covers, the measured effect, that the calls have to be qualified, and that hypot, sqrt and abs are deliberately absent because implementations already evaluate those with SIMD instructions. Spells out the two things a user has to decide on: the accuracy trade (4 ULP where the scalar routines stay below 1, errno and the exception flags unspecified, not bit-wise reproducible against a scalar run) and VIR_DISABLE_SIMD_VECMATH as the way out. Also that the glibc which decides whether anything is routed is the one you build against, not the one you run on, since a toolchain with an old sysroot forwards even on a recent host. On C++26: std::simd specifies these in namespace std::simd and implementations are expected to vectorize them, so this header is unnecessary there. It is scoped to the <experimental/simd> backend and its overloads are constrained to stdx::simd, so it cannot interpose on std::simd; the README says to call std::simd::sinh directly once the standard library provides it. Noted in the header as well, with the condition under which that would have to be revisited. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The functions live in vir::vecmath rather than vir::stdx because a declaration there hides the underlying overloads of that name. Nothing guarded that, and every way it broke before was silent: pow(x, 2.0) stopped compiling, every one-argument name became ambiguous under `using namespace vir::stdx;`, and the two-argument spellings resolved to the underlying overload and quietly took the slow path. The new file asserts the opposite of each: unqualified calls through a using-directive, all the two-argument spellings, hypot in both arities, and that vir::stdx's values are still the underlying implementation's exactly, rather than a vector math library's few ULP away. It needs its own translation unit for the using-directive at namespace scope, which is the check that would otherwise leak into the other two tests. Worth the file: pulling a single function into vir::stdx, the smallest form of the regression, is missed by both existing tests and caught only by this one. Pulling in the whole set is caught by both. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
3 tasks
ax3l
changed the title
simd_vecmath.h: transcendentals via a vector math library
ax3l
commented
Aug 25, 2026
ax3l
commented
Aug 25, 2026
ax3l
commented
Aug 25, 2026
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
The test named exp, cbrt and fabs, which vir's own simd implementation does not provide. Where <experimental/simd> is absent it is that implementation which gets used, so the file failed to compile on gcc 9 to 12, clang 20, every clang-libcxx job and emscripten. This is the same mistake the review caught in the other two files, made again in a file written after it. Use tanh, tan, log2 and erf instead, which both backends have, and drop fabs; abs already covers that check. The test still catches the single-function interposition it exists for. Verified in both configurations this time, which is what would have caught it: -DVIR_DISABLE_STDX_SIMD for vir's own simd, and the default for libstdc++'s. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
vir's own simd, used where the standard library has no <experimental/simd>, implements most of the Parallelism TS math set but not these four. sin, cos, log, log2, erf, pow and hypot are all there; exp is not, which is a surprising hole to land in. It surfaces as soon as anything portable calls them. AMReX builds on macOS, where Apple Clang and libc++ select this simd, failed to compile its amrex::Math SIMD overloads with error: no member named 'exp' in namespace 'vir::stdx' for exp, exp2, expm1 and cbrt, and for nothing else. They go in with the same SIMD_MATH_1ARG the neighbouring functions use, so they are per-lane like the rest of this fallback. Nothing here is about the vector math library: a target that reaches this code has no libmvec to route to in the first place. Worth noting where this does NOT belong. The gap is in the simd type itself, so filling it in vir::vecmath, or in each consumer, would leave vir::stdx incomplete for everyone calling it directly.
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
My naïve attempt to survive with
std::simd+vir-simdin C++20 using math that is more than additions and multiplications #52:Adds
vir/simd_vecmath.h: the transcendentals asvir::vecmath::sinh(x), calling glibc's libmvec entry points directly by their x86-64 vector-ABI names. No compiler flags, no auto-vectorization, inert wherever libmvec is not reachable.Per 2^20 evaluations, width 4, GCC 13.3,
-O3 -march=native, i9-12900H, idle, best of 7:vir::stdxvir::vecmathsinhexpsinsin/cosgain least because libstdc++ already vectorizes those two; everything else it evaluates one lane at a time. In ImpactX an element push calling sin/cos/sinh/cosh per particle got 2.4x faster, the whole run 1.7x. Design, caveats and test coverage: #52Is this the right approach? How would you suggest do I best bridge the math-function-gap until C++26 @mattkretz ? :)