diff --git a/Makefile b/Makefile index 0b06186..4f4504c 100644 --- a/Makefile +++ b/Makefile @@ -4,6 +4,9 @@ # Tests for vir-simd extensions to std::experimental::simd ext_tests = for_each \ + simd_vecmath \ + simd_vecmath_disabled \ + simd_vecmath_no_interposition \ transform \ transform_reduce diff --git a/README.md b/README.md index 8b4b049..e95910f 100644 --- a/README.md +++ b/README.md @@ -32,6 +32,8 @@ gains a proper implementation provides. + [Usable algorithms](#usable-algorithms) + [Example](#example) + [Execution policy modifiers](#execution-policy-modifiers) + - [Vectorized transcendental + math](#vectorized-transcendental-math) - [Bitwise operators for floating-point `simd`](#bitwise-operators-for-floating-point-simd) - [Conversion between `std::bitset` and @@ -289,6 +291,52 @@ its behavior: Determine from run-time information (i.e. add a branch) whether a prologue for alignment of the main chunked iteration might be more efficient. +### Vectorized transcendental math + +```c++ +#include + +auto y = vir::vecmath::sinh(x); +``` + +SIMD hardware has instructions for `sqrt` and `abs` but not for the +transcendental functions, so `simd` implementations evaluate `sin`, `exp`, +`sinh` and friends one lane at a time. A vectorized kernel calling them can +therefore end up slower than the scalar one it replaced. Where glibc's +`libmvec` is available, the functions in `vir::vecmath` hand a whole register +to it instead, by calling its entry points directly. + +Covered are: `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `atan2`, `sinh`, `cosh`, +`tanh`, `asinh`, `acosh`, `atanh`, `exp`, `exp2`, `expm1`, `log`, `log2`, +`log10`, `log1p`, `pow`, `cbrt`, `erf` and `erfc`. + +Call them qualified. An unqualified `sinh(x)` on a `simd` argument resolves to +the underlying implementation through argument-dependent lookup, and no +using-declaration changes that. + +`hypot`, `sqrt` and `abs` are deliberately absent: implementations already +evaluate those with SIMD instructions, so use `vir::stdx` for them. + +Every name is declared in every configuration. Where no vector math library is +reachable — a glibc older than 2.22 (use 2.35+ for full coverage of the list above), +a different libc, another architecture, or +a build with `VIR_DISABLE_SIMD_VECMATH` defined — the call forwards to the +underlying implementation and behaves exactly as before. Note that the relevant +glibc is the one you *build* against, so a toolchain with an old sysroot +forwards even on a recent host. + +> **Accuracy.** Vector math libraries trade accuracy for speed: glibc documents +> a maximum error of 4 ULP for `libmvec`, where its scalar routines stay below +> 1 ULP, and `errno` and the floating-point exception flags become unspecified. +> Results are therefore not bit-wise reproducible against a scalar evaluation. +> Define `VIR_DISABLE_SIMD_VECMATH` if your application needs either. + +> **C++26.** `std::simd` specifies these functions in `namespace std::simd` +> ([\[simd.math\]](https://eel.is/c++draft/simd#math)) and implementations are +> expected to vectorize them, which makes this header unnecessary. It targets +> the `` backend and never interposes on `std::simd`, so +> call `std::simd::sinh` directly once your standard library provides it. + ### Bitwise operators for floating-point `simd` ```c++ diff --git a/testsuite/tests/simd_vecmath.cc b/testsuite/tests/simd_vecmath.cc new file mode 100644 index 0000000..1a2ab24 --- /dev/null +++ b/testsuite/tests/simd_vecmath.cc @@ -0,0 +1,312 @@ +/* SPDX-License-Identifier: GPL-3.0-or-later */ +/* Copyright © 2026 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH + * Matthias Kretz + * Copyright © 2026 Axel Huebl + */ + +// only: float|double|ldouble * * * +// expensive: * [1-9] * * +#include "bits/main.h" +#include + +#include +#include + +/* Coverage for vir/simd_vecmath.h + * + * Everything here is inside VIR_HAVE_SIMD_VECMATH. Where the header is inert + * -- another architecture, another libc, a glibc too old, vir's own simd, or + * the opt-out -- vir::vecmath does not exist and there is nothing of this + * header's to test. Testing the underlying implementation instead would only + * measure how well libstdc++ reduces arguments, which it does badly enough to + * fail (its vector cos returns inf for finite_max). + */ +#ifdef VIR_HAVE_SIMD_VECMATH + +/* Which functions actually reach the vector math library depends on the glibc + * that declared them, so the tolerance has to follow the same split rather + * than the simd's width alone. + */ +#ifdef VIR_HAVE_SIMD_VECMATH_EXTENDED +constexpr bool extended_routed = true; +#else +constexpr bool extended_routed = false; +#endif + +#define VECMATH_TESTER(name_) \ + make_tester("vir::vecmath::" #name_, \ + [](auto... xs) { return vir::vecmath::name_(xs...); }, \ + [](auto... xs) { return std::name_(xs...); }, __FILE__, __LINE__) + +template + int + ulp_distance(T a, T b) + { + constexpr int far = 1 << 30; + if (a == b or (std::isnan(a) and std::isnan(b))) + return 0; + if (std::isnan(a) != std::isnan(b)) + return far; + + static_assert(sizeof(T) == 4 or sizeof(T) == 8, "no integer type to step through"); + using U = std::conditional_t; + constexpr U sign = U(1) << (sizeof(U) * 8 - 1); + U ua, ub; + std::memcpy(&ua, &a, sizeof(T)); + std::memcpy(&ub, &b, sizeof(T)); + + // monotonic unsigned key, so that a plain difference counts steps + const auto key = [sign](U u) { return (u & sign) ? U(~u) : U(u | sign); }; + const U ka = key(ua); + const U kb = key(ub); + const U d = ka > kb ? ka - kb : kb - ka; + return d > U(far) ? far : int(d); + } + +/* Width sweep + * + * Which path a simd takes, and how many chunks the vector math library is + * called with, is decided entirely by its width. The harness instantiates only + * the scalar and the native ABI unless the expensive tests are enabled, and + * its ULP helper does not work with fixed_size ABIs on top of libstdc++'s + * simd, so the widths are swept here, lane by lane, with nothing from the + * harness involved. + */ +template + void + test_one_width(const char* name, bool routed, FSimd&& fsimd, FScalar&& fscalar, + std::initializer_list inputs) + { + using V = vir::stdx::fixed_size_simd; + /* 4 ULP is what libmvec documents. Where the call is not routed the answer + * comes from the underlying implementation, which is not the scalar + * routine either: libstdc++ carries its own sin and cos, good to 1 ULP. + */ + const int allowed + = (routed and vir::vecmath_detail::use_vecmath) ? 4 : 1; + + alignas(vir::stdx::memory_alignment_v) T lane[Width]; + auto it = inputs.begin(); + for (int i = 0; i < Width; ++i, ++it) + { + if (it == inputs.end()) + it = inputs.begin(); + lane[i] = *it; + } + + V x; + x.copy_from(lane, vir::stdx::vector_aligned); + const V got = fsimd(x); + + for (int i = 0; i < Width; ++i) + { + const T expect = fscalar(lane[i]); + const int d = ulp_distance(T(got[i]), expect); + VERIFY(d <= allowed) + << name << " at width " << Width << ", lane " << i << ": got " << T(got[i]) + << ", expected " << expect << " (" << d << " ULP, allowed " << allowed << ')'; + } + } + +/* How wide a chunk is decides how many calls into the vector math library a + * simd costs. Narrowing it keeps every result correct, so no value comparison + * can notice; assert the selection directly instead. + */ +template + void + check_chunk_width() + { + using vir::vecmath_detail::chunk_width; + using vir::vecmath_detail::is_supported_width; + using vir::vecmath_detail::native_lanes; + + if constexpr (is_supported_width) + { + constexpr int chunk = chunk_width::value; + constexpr int native = native_lanes::value; + constexpr int smallest = 16 / int(sizeof(T)); + + static_assert(Width % chunk == 0, + "a chunk that does not divide the width would run off the end"); + static_assert(chunk <= native, + "a chunk wider than the target's registers is not callable"); + static_assert(chunk >= smallest, + "the narrowest entry point takes a full 128-bit register"); + static_assert(Width % native != 0 or chunk == native, + "a width the native register divides has to use the native chunk"); + } + } + +#define SWEEP_1(name_, routed_, values_) \ + (test_one_width(#name_, routed_, [](auto v) { return vir::vecmath::name_(v); }, \ + [](T v) { return std::name_(v); }, values_), ...) + +#define SWEEP_2(name_, routed_, second_, values_) \ + (test_one_width(#name_, routed_, \ + [](auto v) { return vir::vecmath::name_(v, T(second_)); }, \ + [](T v) { return std::name_(v, T(second_)); }, values_), ...) + +template + void + sweep_widths() + { + (check_chunk_width(), ...); + + // a spread of values, cycled to fill each width, inside the domain of each group + const std::initializer_list general + = {T(0.5), T(-0.5), T(1), T(-1), T(0), T(2), T(-2), T(0.25)}; + const std::initializer_list positive + = {T(0.5), T(1), T(1.5), T(2), T(3), T(0.25), T(10), T(1.25)}; + const std::initializer_list unit + = {T(0.5), T(-0.5), T(1), T(-1), T(0), T(0.25), T(-0.25), T(0.75)}; + const std::initializer_list above_one + = {T(1), T(1.5), T(2), T(3), T(10), T(1.25), T(5), T(100)}; + + constexpr bool base = true; // sin cos exp log pow, since glibc 2.22 + constexpr bool ext = extended_routed; // the rest, since glibc 2.35 + + SWEEP_1(sin, base, general); SWEEP_1(cos, base, general); + SWEEP_1(exp, base, general); SWEEP_1(log, base, positive); + SWEEP_2(pow, base, 2.5, positive); + + SWEEP_1(tan, ext, general); SWEEP_1(atan, ext, general); + SWEEP_1(sinh, ext, general); SWEEP_1(cosh, ext, general); + SWEEP_1(tanh, ext, general); SWEEP_1(asinh, ext, general); + SWEEP_1(cbrt, ext, general); SWEEP_1(erf, ext, general); + SWEEP_1(erfc, ext, general); SWEEP_1(exp2, ext, general); + SWEEP_1(expm1, ext, general); SWEEP_1(asin, ext, unit); + SWEEP_1(acos, ext, unit); SWEEP_1(atanh, ext, unit); + SWEEP_1(acosh, ext, above_one); SWEEP_1(log2, ext, positive); + SWEEP_1(log10, ext, positive); SWEEP_1(log1p, ext, positive); + SWEEP_2(atan2, ext, 2, general); + } + +#undef SWEEP_1 +#undef SWEEP_2 +#endif // VIR_HAVE_SIMD_VECMATH + +template + void + test() + { +#ifdef VIR_HAVE_SIMD_VECMATH + using T = typename V::value_type; + using Abi = typename V::abi_type; + + /* test_values reaches vir::detail::bit_cast through the harness's ULP + * helper, which does not accept fixed_size ABIs on top of libstdc++'s + * simd. Those widths are covered by the sweep instead. + */ + constexpr bool harness_usable + = !std::is_same_v>; + + constexpr bool vecmath = vir::vecmath_detail::use_vecmath; + + if constexpr (harness_usable) + { + constexpr T inf = vir::infinity_v; + constexpr T nan = vir::quiet_NaN_v; + constexpr T denorm_min = vir::denorm_min_v; + constexpr T norm_min = vir::norm_min_v; + constexpr T max = vir::finite_max_v; + + const std::initializer_list edge_values + = {+0., -0., 1., -1., 0.5, -0.5, 2., -2., inf, -inf, nan, + denorm_min, -denorm_min, norm_min, norm_min / 3, max, -max}; + + // libmvec does not reproduce the scalar routines' exception flags + FloatExceptCompare::ignore = vecmath; + + // the functions glibc has had since 2.22 + vir::test::setFuzzyness(vecmath ? 4 : 1); + vir::test::setFuzzyness(vecmath ? 4 : 1); + vir::test::setFuzzyness(1); + + test_values(edge_values, {5000}, + VECMATH_TESTER(sin), VECMATH_TESTER(cos), VECMATH_TESTER(exp)); + test_values({norm_min, denorm_min, 0.5, 1., 2., max, inf, nan}, {5000, denorm_min, max}, + VECMATH_TESTER(log)); + + // and the ones it grew in 2.35 + const bool ext = vecmath and extended_routed; + FloatExceptCompare::ignore = ext; + vir::test::setFuzzyness(ext ? 4 : 1); + vir::test::setFuzzyness(ext ? 4 : 1); + + test_values(edge_values, {5000}, + VECMATH_TESTER(tan), VECMATH_TESTER(atan), + VECMATH_TESTER(erf), VECMATH_TESTER(erfc)); + test_values(edge_values, {5000}, + VECMATH_TESTER(sinh), VECMATH_TESTER(cosh), VECMATH_TESTER(tanh), + VECMATH_TESTER(asinh), VECMATH_TESTER(cbrt)); + test_values(edge_values, {5000}, + VECMATH_TESTER(exp2), VECMATH_TESTER(expm1)); + test_values({-1., -0.5, +0., -0., 0.5, 1., denorm_min, nan}, {5000, T(-1), T(1)}, + VECMATH_TESTER(asin), VECMATH_TESTER(acos), VECMATH_TESTER(atanh)); + test_values({1., 1.5, 2., max, inf, nan}, {5000, T(1), max}, + VECMATH_TESTER(acosh)); + test_values({norm_min, denorm_min, 0.5, 1., 2., max, inf, nan}, {5000, denorm_min, max}, + VECMATH_TESTER(log2), VECMATH_TESTER(log10)); + test_values({-1., -0.5, +0., -0., 0.5, 1., max, inf, nan}, {5000, T(-1), max}, + VECMATH_TESTER(log1p)); + test_values_2arg({+0., -0., 0.5, 1., 2., 3., inf, -inf, nan, norm_min, max}, + {5000}, VECMATH_TESTER(atan2)); + + FloatExceptCompare::ignore = vecmath; + vir::test::setFuzzyness(vecmath ? 4 : 1); + vir::test::setFuzzyness(vecmath ? 4 : 1); + test_values_2arg({+0., -0., 0.5, 1., 2., 3., -2., inf, -inf, nan, norm_min}, + {2000, T(0), T(10)}, VECMATH_TESTER(pow)); + + FloatExceptCompare::ignore = false; + vir::test::setFuzzyness(0); + vir::test::setFuzzyness(0); + } + + // every width that decides chunking, once per element type + if constexpr (std::is_same_v + and (sizeof(T) == 4 or sizeof(T) == 8)) + sweep_widths(); + + /* The two-argument overloads must accept a scalar on either side. Deducing + * both parameters instead of only the first compiles fine and silently + * takes these spellings away from every caller. The exponent is not a + * whole number on purpose: with 2 the vector and scalar results agree bit + * for bit, which makes the comparison vacuous. + */ + { + const V x = make_value_unknown(V([](auto i) { return T(1) + T(i) * T(0.25); })); + const V e = make_value_unknown(V(T(2.5))); + const T e_scalar = make_value_unknown(T(2.5)); + + COMPARE(vir::vecmath::pow(x, e_scalar), vir::vecmath::pow(x, e)); + COMPARE(vir::vecmath::pow(e_scalar, x), vir::vecmath::pow(e, x)); + COMPARE(vir::vecmath::atan2(x, e_scalar), vir::vecmath::atan2(x, e)); + COMPARE(vir::vecmath::atan2(e_scalar, x), vir::vecmath::atan2(e, x)); + + // an argument that merely converts to the element type has to work too + COMPARE(vir::vecmath::pow(x, make_value_unknown(2)), vir::vecmath::pow(x, V(T(2)))); + COMPARE(vir::vecmath::pow(make_value_unknown(2), x), vir::vecmath::pow(V(T(2)), x)); + } + + /* These live in their own namespace now, so vir::stdx keeps every overload + * it had. That is the whole point of not declaring them there. + */ + { + COMPARE(vir::stdx::hypot(V(T(3)), V(T(4))), V(T(5))); + VERIFY(all_of(vir::stdx::hypot(V(T(3)), V(T(4)), V(T(0))) == V(T(5)))); + COMPARE(vir::stdx::sqrt(V(T(4))), V(T(2))); + COMPARE(vir::stdx::abs(V(T(-2))), V(T(2))); + COMPARE(vir::stdx::fabs(V(T(-2))), V(T(2))); + COMPARE(vir::stdx::pow(V(T(2)), T(3)), V(T(8))); // the scalar-broadcast form + COMPARE(vir::stdx::pow(T(2), V(T(3))), V(T(8))); // and the reversed one + } + + // taking an address forces the out-of-line copy the ISA tag keeps apart + { + using Fn = V (*)(const V&); + const Fn f = &vir::vecmath::sin; + COMPARE(f(V(T(0))), V(T(0))); + } +#endif // VIR_HAVE_SIMD_VECMATH + } diff --git a/testsuite/tests/simd_vecmath_disabled.cc b/testsuite/tests/simd_vecmath_disabled.cc new file mode 100644 index 0000000..dbc6b9c --- /dev/null +++ b/testsuite/tests/simd_vecmath_disabled.cc @@ -0,0 +1,56 @@ +/* SPDX-License-Identifier: GPL-3.0-or-later */ +/* Copyright © 2026 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH + * Matthias Kretz + * Copyright © 2026 Axel Huebl + */ + +// only: float|double|ldouble * * * +// expensive: * [1-9] * * + +/* VIR_DISABLE_SIMD_VECMATH has to switch the header off completely + * + * A user who needs errno, the exception flags or sub-ULP accuracy has to be + * able to opt out even when the header arrives through another one. The + * compile-time check below is the substance of it; the rest confirms that + * vir::stdx is left exactly as it was, which is also what the header claims + * when it is enabled, since it no longer declares anything there. + */ +#define VIR_DISABLE_SIMD_VECMATH 1 +#include "bits/main.h" +#include + +#ifdef VIR_HAVE_SIMD_VECMATH +#error "VIR_DISABLE_SIMD_VECMATH did not disable vir/simd_vecmath.h" +#endif + +template + void + test() + { + using T = typename V::value_type; + + /* Only names the underlying implementation is known to provide: with vir's + * own simd, exp, exp2, expm1, cbrt and fabs do not exist at all, and this + * file is compiled in that configuration too. + */ + const V x = make_value_unknown(V([](auto i) { return T(0.25) + T(i) * T(0.125); })); + const V unit = make_value_unknown(V([](auto i) { return T(i) / T(V::size()); })); + + // still reachable, and still answering correctly + COMPARE(vir::stdx::sin(V(T(0))), V(T(0))); + COMPARE(vir::stdx::cos(V(T(0))), V(T(1))); + COMPARE(vir::stdx::sinh(V(T(0))), V(T(0))); + COMPARE(vir::stdx::log(V(T(1))), V(T(0))); + COMPARE(vir::stdx::asin(V(T(0))), V(T(0))); + VERIFY(all_of(vir::stdx::sin(x) * vir::stdx::sin(x) + + vir::stdx::cos(x) * vir::stdx::cos(x) > V(T(0.99)))); + VERIFY(all_of(vir::stdx::asin(unit) <= V(T(1.5708)))); + + // and the full overload set is intact + COMPARE(vir::stdx::pow(V(T(2)), V(T(3))), V(T(8))); + COMPARE(vir::stdx::pow(V(T(2)), T(3)), V(T(8))); + COMPARE(vir::stdx::pow(T(2), V(T(3))), V(T(8))); + COMPARE(vir::stdx::hypot(V(T(3)), V(T(4))), V(T(5))); + COMPARE(vir::stdx::sqrt(V(T(4))), V(T(2))); + COMPARE(vir::stdx::abs(V(T(-2))), V(T(2))); + } diff --git a/testsuite/tests/simd_vecmath_no_interposition.cc b/testsuite/tests/simd_vecmath_no_interposition.cc new file mode 100644 index 0000000..3ae194c --- /dev/null +++ b/testsuite/tests/simd_vecmath_no_interposition.cc @@ -0,0 +1,89 @@ +/* SPDX-License-Identifier: GPL-3.0-or-later */ +/* Copyright © 2026 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH + * Matthias Kretz + * Copyright © 2026 Axel Huebl + */ + +// only: float|double|ldouble * * * +// expensive: * [1-9] * * + +/* Including vir/simd_vecmath.h must leave vir::stdx exactly as it was + * + * The functions could have been declared in vir::stdx, which would have made + * vir::stdx::sinh resolve to them. They are not, because a declaration there + * hides the underlying overloads of that name, and everything below is what + * that would have cost. Each check fails to *compile* rather than to run if + * the header ever starts interposing again, which is the point: every one of + * these was silent when it broke. + * + * This needs its own translation unit for the using-directive at namespace + * scope, which is the check that would otherwise leak into the other tests. + */ +#include "bits/main.h" +#include + +// the directive that turns an interposing overload set into an ambiguity +using namespace vir::stdx; + +template + void + test() + { + using T = typename V::value_type; + + const V x = make_value_unknown(V([](auto i) { return T(1) + T(i) * T(0.125); })); + const V two = make_value_unknown(V(T(2))); + const T two_scalar = make_value_unknown(T(2)); + + /* Unqualified, through the using-directive above. If vir::vecmath's + * overloads were visible here too, every one of these would be ambiguous. + */ + /* Only names both backends provide: vir's own simd implementation, used + * where is absent, has no exp, exp2, expm1, cbrt or + * fabs, and this file is compiled in that configuration too. + */ + VERIFY(all_of(sinh(x) == vir::stdx::sinh(x))) << "unqualified sinh"; + VERIFY(all_of(cosh(x) == vir::stdx::cosh(x))) << "unqualified cosh"; + VERIFY(all_of(tanh(x) == vir::stdx::tanh(x))) << "unqualified tanh"; + VERIFY(all_of(sin(x) == vir::stdx::sin(x))) << "unqualified sin"; + VERIFY(all_of(cos(x) == vir::stdx::cos(x))) << "unqualified cos"; + VERIFY(all_of(tan(x) == vir::stdx::tan(x))) << "unqualified tan"; + VERIFY(all_of(log(x) == vir::stdx::log(x))) << "unqualified log"; + VERIFY(all_of(log2(x) == vir::stdx::log2(x))) << "unqualified log2"; + VERIFY(all_of(atan(x) == vir::stdx::atan(x))) << "unqualified atan"; + VERIFY(all_of(erf(x) == vir::stdx::erf(x))) << "unqualified erf"; + + /* The two-argument overload set. libstdc++ generates two overloads per + * function, the first with a second parameter excluded from deduction, + * which is what lets a scalar be broadcast. Reproducing that set wrongly + * is how pow(x, 2.0) stopped compiling once before. + */ + VERIFY(all_of(pow(x, two) == vir::stdx::pow(x, two))) << "pow(simd, simd)"; + VERIFY(all_of(pow(x, two_scalar) == vir::stdx::pow(x, two))) << "pow(simd, scalar)"; + VERIFY(all_of(pow(two_scalar, x) == vir::stdx::pow(two, x))) << "pow(scalar, simd)"; + VERIFY(all_of(atan2(x, two_scalar) == vir::stdx::atan2(x, two))) << "atan2(simd, scalar)"; + VERIFY(all_of(atan2(two_scalar, x) == vir::stdx::atan2(two, x))) << "atan2(scalar, simd)"; + + // an argument that merely converts to the element type + VERIFY(all_of(pow(x, make_value_unknown(2)) == vir::stdx::pow(x, two))) << "pow(simd, int)"; + + /* Names the header deliberately leaves alone. hypot in particular carries + * a three-argument form and converting overloads that a declaration in + * vir::stdx would have taken away. + */ + COMPARE(hypot(V(T(3)), V(T(4))), V(T(5))); + VERIFY(all_of(hypot(V(T(3)), V(T(4)), V(T(0))) == V(T(5)))) << "hypot, three arguments"; + COMPARE(sqrt(V(T(4))), V(T(2))); + COMPARE(abs(V(T(-2))), V(T(2))); + + /* And the values are the underlying implementation's, not a vector math + * library's: identical, not merely close. Anything routed through + * vir::vecmath would differ here by the few ULP that costs. + */ + for (std::size_t i = 0; i < V::size(); ++i) + { + const T xi = T(x[i]); + COMPARE(T(vir::stdx::sinh(x)[i]), std::sinh(xi)) << "lane " << i; + COMPARE(T(vir::stdx::tanh(x)[i]), std::tanh(xi)) << "lane " << i; + } + } diff --git a/vir/simd.h b/vir/simd.h index 5f2d6a8..c77d144 100644 --- a/vir/simd.h +++ b/vir/simd.h @@ -2595,6 +2595,12 @@ namespace vir::stdx SIMD_MATH_1ARG(log2, simd) SIMD_MATH_1ARG(logb, simd) + // exponentials and roots + SIMD_MATH_1ARG(exp, simd) + SIMD_MATH_1ARG(exp2, simd) + SIMD_MATH_1ARG(expm1, simd) + SIMD_MATH_1ARG(cbrt, simd) + #undef SIMD_MATH_1ARG #undef SIMD_MATH_1ARG_FIXED #undef SIMD_MATH_2ARG diff --git a/vir/simd_vecmath.h b/vir/simd_vecmath.h new file mode 100644 index 0000000..f28ab6c --- /dev/null +++ b/vir/simd_vecmath.h @@ -0,0 +1,579 @@ +/* SPDX-License-Identifier: LGPL-3.0-or-later */ +/* Copyright © 2026 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH + * Matthias Kretz + * Copyright © 2026 Axel Huebl + */ + +#ifndef VIR_SIMD_VECMATH_H_ +#define VIR_SIMD_VECMATH_H_ + +#include "simd.h" + +#include +#include +#include + +/* Transcendental math functions evaluated by a vector math library. + * + * SIMD hardware has instructions for sqrt and abs, but not for the + * transcendental functions. std::experimental::simd therefore falls back to + * calling the scalar libm function once per lane, which makes a vectorized + * kernel calling sin, exp or sinh slower than the scalar one it replaced. + * + * A vector math library computes a whole register worth of results per call. + * This header offers the functions glibc's libmvec provides as + * vir::vecmath::sinh(x) and so on, and leaves everything else to the + * underlying simd implementation. Call them qualified: an unqualified sinh(x) + * resolves to the underlying implementation through argument-dependent + * lookup, and no using-declaration changes that. + * + * The calls go to the vector entry points directly, using the x86-64 vector + * function ABI names, so no compiler flags and no auto-vectorization are + * involved. + * + * What this costs. 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 here: + * + * - errno may be set where the scalar routine leaves it alone, and left alone + * where the scalar routine sets it, + * - the exception flags may gain FE_INEXACT on exact results, may miss flags + * the scalar routine raises, and may raise a different one (atanh(1) reports + * FE_INVALID rather than only FE_DIVBYZERO), + * - flags and errno are per call, not per lane, so one lane hitting a pole + * leaves them set for the whole vector, + * - results are accurate to 4 ULP where the scalar routines stay below 1, so + * they are not bit-wise identical to a scalar evaluation. + * + * Sign of zero, infinities, NaNs and denormals are handled the same as by the + * scalar routines. Define VIR_DISABLE_SIMD_VECMATH to keep the underlying + * implementation, which has none of the above caveats. + * + * On C++26. std::simd specifies these functions in namespace std::simd + * ([simd.math]) and implementations are expected to vectorize them, which makes + * this header unnecessary there. It is scoped to the + * backend by the VIR_HAVE_STD_SIMD condition above, and its overloads are + * constrained to stdx::simd, so it never interposes on std::simd. Should + * vir::stdx ever be backed by std::simd, this header has to be revisited: the + * routing below would then be shadowing the standard library's own vectorized + * math rather than filling a gap in it. + * + * For what the standard intends here, see P1928R15 section 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 reflected in the wording. Dropping errno is therefore the direction of + * travel; the exception flags are where a vector math library falls short of + * it. No accuracy bound is specified either way, so the ULP figures above are + * glibc's own documentation rather than anything guaranteed. + */ + +// __GLIBC__ only exists once a libc header has been seen, so pull it in first +#if __has_include() +#include +#endif + +// __ILP32__ excludes x32, where __x86_64__ is defined but the ABI is not this one +#if defined __x86_64__ && !defined __ILP32__ && defined __GLIBC__ \ + && defined VIR_HAVE_STD_SIMD && !defined VIR_DISABLE_SIMD_VECMATH +#ifdef __GLIBC_PREREQ + +#if __GLIBC_PREREQ(2, 22) +#define VIR_HAVE_SIMD_VECMATH 1 +#endif +#if __GLIBC_PREREQ(2, 35) +//! glibc 2.35 grew vector variants for everything beyond sin, cos, exp, log and pow +#define VIR_HAVE_SIMD_VECMATH_EXTENDED 1 +#endif + +#endif // __GLIBC_PREREQ +#endif // __x86_64__ && __GLIBC__ && VIR_HAVE_STD_SIMD + +#ifdef VIR_HAVE_SIMD_VECMATH + +namespace vir::vecmath_detail +{ + template + using vec [[gnu::vector_size(Width * sizeof(T))]] = T; + + using v2d = vec; + using v4d = vec; + using v8d = vec; + using v4f = vec; + using v8f = vec; + using v16f = vec; + + /* Which lane counts can be handed to the vector math library + * + * The ISA letter in the symbol name says which instruction set the callee + * uses: b is SSE2, c is AVX, d is AVX2 and e is AVX-512. A 256-bit call is + * only available as AVX2 when the translation unit is built for AVX2, so the + * letter follows the compiled-for ISA, while the lane count follows the + * register width being passed. + */ + template + struct native_lanes + { +#if defined __AVX512F__ + static constexpr int value = 64 / int(sizeof(T)); +#elif defined __AVX__ + static constexpr int value = 32 / int(sizeof(T)); +#else + static constexpr int value = 16 / int(sizeof(T)); +#endif + }; + + //! true if a simd of Width elements of T can be evaluated in whole chunks + template + inline constexpr bool is_supported_width + = (std::is_same_v || std::is_same_v) + && Width >= (16 / int(sizeof(T))) + && Width % (16 / int(sizeof(T))) == 0; + + /* Chunk width used for a simd of Width elements + * + * The largest chunk the ISA supports that divides Width, so that a simd + * wider than one register is evaluated in a few full calls rather than + * falling back to scalar. + */ + template + struct chunk_width + { + static constexpr int native = native_lanes::value; + static constexpr int value + = Width % native == 0 ? native + : (native > 16 / int(sizeof(T)) && Width % (native / 2) == 0) ? native / 2 + : 16 / int(sizeof(T)); + }; +} // namespace vir::vecmath_detail + +/* Naming the libmvec entry points + * + * The x86-64 vector function ABI spells them _ZGV_, + * where the isa letter says which instruction set the callee uses: b is SSE2, + * c is AVX, d is AVX2, e is AVX-512. A 256-bit call therefore has two spellings + * and the right one is the one matching what this translation unit is built + * for, while the lane count follows the register being passed. + * + * The c (AVX) entry points are the one class glibc does not resolve through an + * ifunc: they are wrappers that split the argument and call the SSE2 routine + * twice. They are still far ahead of one scalar call per lane, but a target + * without AVX2 should not expect a true 256-bit routine. + */ +#if defined __AVX2__ +# define VIR_VECMATH_ISA256 d +#else +# define VIR_VECMATH_ISA256 c +#endif + +/* The chunk callers below are always_inline, so normally no out-of-line copy + * exists, but -fkeep-inline-functions or taking an address emits one, and its + * body differs per instruction set exactly as the public overloads' does. An + * inline namespace named after the ISA keeps those symbols apart without + * touching how they are called. + */ +#if defined __AVX512F__ +# define VIR_VECMATH_ISA_NS isa_avx512 +#elif defined __AVX2__ +# define VIR_VECMATH_ISA_NS isa_avx2 +#elif defined __AVX__ +# define VIR_VECMATH_ISA_NS isa_avx +#else +# define VIR_VECMATH_ISA_NS isa_sse2 +#endif + +#define VIR_VECMATH_CAT_(a, b) a##b +#define VIR_VECMATH_CAT(a, b) VIR_VECMATH_CAT_(a, b) +#define VIR_VECMATH_SYM(isa, rest) VIR_VECMATH_CAT(VIR_VECMATH_CAT(_ZGV, isa), rest) + +/* Declaring them + * + * Only the widths the target ISA actually has: declaring a function that + * returns a 512-bit vector without AVX-512 enabled would change the ABI, which + * GCC rightly warns about (-Wpsabi). + */ +#if defined __AVX512F__ +# define VIR_VECMATH_DECL_1_512(name) \ + vir::vecmath_detail::v8d _ZGVeN8v_##name (vir::vecmath_detail::v8d); \ + vir::vecmath_detail::v16f _ZGVeN16v_##name##f (vir::vecmath_detail::v16f); +# define VIR_VECMATH_DECL_2_512(name) \ + vir::vecmath_detail::v8d _ZGVeN8vv_##name (vir::vecmath_detail::v8d, \ + vir::vecmath_detail::v8d); \ + vir::vecmath_detail::v16f _ZGVeN16vv_##name##f (vir::vecmath_detail::v16f, \ + vir::vecmath_detail::v16f); +#else +# define VIR_VECMATH_DECL_1_512(name) +# define VIR_VECMATH_DECL_2_512(name) +#endif + +#if defined __AVX__ +# define VIR_VECMATH_DECL_1_256(name) \ + vir::vecmath_detail::v4d VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N4v_##name) \ + (vir::vecmath_detail::v4d); \ + vir::vecmath_detail::v8f VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N8v_##name##f) \ + (vir::vecmath_detail::v8f); +# define VIR_VECMATH_DECL_2_256(name) \ + vir::vecmath_detail::v4d VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N4vv_##name) \ + (vir::vecmath_detail::v4d, vir::vecmath_detail::v4d); \ + vir::vecmath_detail::v8f VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N8vv_##name##f) \ + (vir::vecmath_detail::v8f, vir::vecmath_detail::v8f); +#else +# define VIR_VECMATH_DECL_1_256(name) +# define VIR_VECMATH_DECL_2_256(name) +#endif + +#define VIR_VECMATH_DECL_1(name) \ + extern "C" { \ + vir::vecmath_detail::v2d _ZGVbN2v_##name (vir::vecmath_detail::v2d); \ + vir::vecmath_detail::v4f _ZGVbN4v_##name##f (vir::vecmath_detail::v4f); \ + VIR_VECMATH_DECL_1_256(name) \ + VIR_VECMATH_DECL_1_512(name) \ + } + +#define VIR_VECMATH_DECL_2(name) \ + extern "C" { \ + vir::vecmath_detail::v2d _ZGVbN2vv_##name (vir::vecmath_detail::v2d, \ + vir::vecmath_detail::v2d); \ + vir::vecmath_detail::v4f _ZGVbN4vv_##name##f (vir::vecmath_detail::v4f, \ + vir::vecmath_detail::v4f); \ + VIR_VECMATH_DECL_2_256(name) \ + VIR_VECMATH_DECL_2_512(name) \ + } + +/* Selecting the entry point for one chunk + * + * VIR_VECMATH_CALL_1(sin) defines vir::vecmath_detail::call_sin, overloaded on + * the raw vector type, so the chunk loop below simply calls it. + */ +#if defined __AVX512F__ +# define VIR_VECMATH_CALL_1_512(name) \ + VIR_ALWAYS_INLINE v8d call_##name (v8d x) { return _ZGVeN8v_##name(x); } \ + VIR_ALWAYS_INLINE v16f call_##name (v16f x) { return _ZGVeN16v_##name##f(x); } +# define VIR_VECMATH_CALL_2_512(name) \ + VIR_ALWAYS_INLINE v8d call_##name (v8d x, v8d y) { return _ZGVeN8vv_##name(x, y); } \ + VIR_ALWAYS_INLINE v16f call_##name (v16f x, v16f y) \ + { return _ZGVeN16vv_##name##f(x, y); } +#else +# define VIR_VECMATH_CALL_1_512(name) +# define VIR_VECMATH_CALL_2_512(name) +#endif + +#if defined __AVX__ +# define VIR_VECMATH_CALL_1_256(name) \ + VIR_ALWAYS_INLINE v4d call_##name (v4d x) \ + { return VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N4v_##name)(x); } \ + VIR_ALWAYS_INLINE v8f call_##name (v8f x) \ + { return VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N8v_##name##f)(x); } +# define VIR_VECMATH_CALL_2_256(name) \ + VIR_ALWAYS_INLINE v4d call_##name (v4d x, v4d y) \ + { return VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N4vv_##name)(x, y); } \ + VIR_ALWAYS_INLINE v8f call_##name (v8f x, v8f y) \ + { return VIR_VECMATH_SYM(VIR_VECMATH_ISA256, N8vv_##name##f)(x, y); } +#else +# define VIR_VECMATH_CALL_1_256(name) +# define VIR_VECMATH_CALL_2_256(name) +#endif + +/* The chunk callers carry the tag as well: they are always_inline too, but + * -fkeep-inline-functions or taking an address emits them, and the body + * differs per ISA exactly as the public overloads' does. + */ +#define VIR_VECMATH_CALL_1(name) \ + namespace vir::vecmath_detail { inline namespace VIR_VECMATH_ISA_NS { \ + VIR_ALWAYS_INLINE v2d call_##name (v2d x) { return _ZGVbN2v_##name(x); } \ + VIR_ALWAYS_INLINE v4f call_##name (v4f x) { return _ZGVbN4v_##name##f(x); } \ + VIR_VECMATH_CALL_1_256(name) \ + VIR_VECMATH_CALL_1_512(name) \ + } } + +#define VIR_VECMATH_CALL_2(name) \ + namespace vir::vecmath_detail { inline namespace VIR_VECMATH_ISA_NS { \ + VIR_ALWAYS_INLINE v2d call_##name (v2d x, v2d y) { return _ZGVbN2vv_##name(x, y); } \ + VIR_ALWAYS_INLINE v4f call_##name (v4f x, v4f y) { return _ZGVbN4vv_##name##f(x, y); } \ + VIR_VECMATH_CALL_2_256(name) \ + VIR_VECMATH_CALL_2_512(name) \ + } } + +namespace vir::vecmath_detail +{ + /* Discriminator for the emitted symbol + * + * These functions are always_inline, so at -O1 and above no out-of-line copy + * is emitted at all. At -O0, or when the address of one is taken, a weak + * symbol appears, and two translation units built for different instruction + * sets would otherwise agree on its name while disagreeing on its body: the + * linker keeps one, and the other TU ends up calling a libmvec entry point + * its target may not be able to execute. Naming the ISA in the signature + * keeps those symbols apart. libstdc++ solves the same problem the same way, + * see __odr_helper in . + */ + template + struct isa_tag {}; + +#if defined __AVX512F__ + using odr_tag = isa_tag<3>; +#elif defined __AVX2__ + using odr_tag = isa_tag<2>; +#elif defined __AVX__ + using odr_tag = isa_tag<1>; +#else + using odr_tag = isa_tag<0>; +#endif + + /* Keeps a parameter out of template argument deduction + * + * The second argument of a two-argument function must not take part in + * deduction, so that pow(x, 2.0) converts the scalar to a simd instead of + * failing to deduce. This mirrors _Extra_argument_type in libstdc++. + */ + template + struct nondeduced { using type = T; }; + + template + using nondeduced_t = typename nondeduced::type; + + //! true for a first argument that is not a simd but converts to one + template + inline constexpr bool is_convertible_first + = std::is_floating_point_v + && !std::is_same_v, stdx::simd> + && std::is_convertible_v>; + + //! true if the math functions below hand this simd to the vector math library + template + inline constexpr bool use_vecmath + = std::is_floating_point_v + && is_supported_width::size())>; + + //! true if they leave it to the underlying simd implementation instead + template + inline constexpr bool use_fallback + = std::is_floating_point_v + && !is_supported_width::size())>; + + /* Evaluate call on every chunk of x + * + * The round trip through lane[] is what lets this use nothing but the public + * simd interface. It costs nothing: the buffer is a local of exactly the + * chunk alignment, so the stores and loads fold into register moves and the + * generated code is a plain call per chunk. + */ + template + VIR_ALWAYS_INLINE stdx::simd + apply (const stdx::simd& x, F&& call) + { + using V = stdx::simd; + constexpr int width = int(V::size()); + constexpr int chunk = chunk_width::value; + using chunk_type = vec; + + alignas(stdx::memory_alignment_v) T lane[width]; + x.copy_to(lane, stdx::vector_aligned); + + for (int i = 0; i < width; i += chunk) + { + chunk_type v; + std::memcpy(&v, lane + i, sizeof(chunk_type)); + v = call(v); + std::memcpy(lane + i, &v, sizeof(chunk_type)); + } + + V r; + r.copy_from(lane, stdx::vector_aligned); + return r; + } + + //! @see apply + template + VIR_ALWAYS_INLINE stdx::simd + apply (const stdx::simd& x, const stdx::simd& y, F&& call) + { + using V = stdx::simd; + constexpr int width = int(V::size()); + constexpr int chunk = chunk_width::value; + using chunk_type = vec; + + alignas(stdx::memory_alignment_v) T lane_x[width]; + alignas(stdx::memory_alignment_v) T lane_y[width]; + x.copy_to(lane_x, stdx::vector_aligned); + y.copy_to(lane_y, stdx::vector_aligned); + + for (int i = 0; i < width; i += chunk) + { + chunk_type vx, vy; + std::memcpy(&vx, lane_x + i, sizeof(chunk_type)); + std::memcpy(&vy, lane_y + i, sizeof(chunk_type)); + vx = call(vx, vy); + std::memcpy(lane_x + i, &vx, sizeof(chunk_type)); + } + + V r; + r.copy_from(lane_x, stdx::vector_aligned); + return r; + } +} // namespace vir::vecmath_detail + +/* Defining the overloads + * + * The overloads go into vir::stdx, which already pulls in the underlying + * implementation with a using-directive. Qualified lookup stops as soon as it + * finds a declaration in vir::stdx itself, so vir::stdx::sinh means these and + * never the underlying one, without any ambiguity. + * + * That also means the pair below has to cover every simd the underlying + * implementation covers: the fallback overload is not an optimization, it is + * what keeps unsupported element types and widths working. + */ +#define VIR_VECMATH_FN_1(name) \ + VIR_VECMATH_DECL_1(name) \ + VIR_VECMATH_CALL_1(name) \ + namespace vir::vecmath { \ + template \ + VIR_ALWAYS_INLINE \ + std::enable_if_t, stdx::simd> \ + name (const stdx::simd& x) \ + { \ + return vecmath_detail::apply( \ + x, [](auto v) { return vecmath_detail::call_##name(v); }); \ + } \ + \ + template \ + VIR_ALWAYS_INLINE \ + std::enable_if_t, stdx::simd> \ + name (const stdx::simd& x) \ + { return std::experimental::parallelism_v2::name(x); } \ + } + +/* Two-argument functions + * + * The shape follows libstdc++'s _GLIBCXX_SIMD_MATH_CALL2_ exactly, because + * these overloads hide it: a first form whose second parameter is excluded + * from deduction, so that pow(x, 2.0) broadcasts the scalar, and a reversed + * form for pow(2.0, x). Deducing both parameters instead would compile, and + * would silently take those two spellings away from every caller. + */ +#define VIR_VECMATH_FN_2(name) \ + VIR_VECMATH_DECL_2(name) \ + VIR_VECMATH_CALL_2(name) \ + namespace vir::vecmath { \ + template \ + VIR_ALWAYS_INLINE \ + std::enable_if_t, stdx::simd> \ + name (const stdx::simd& x, \ + const vecmath_detail::nondeduced_t>& y) \ + { \ + return vecmath_detail::apply( \ + x, y, [](auto a, auto b) { return vecmath_detail::call_##name(a, b); }); \ + } \ + \ + template \ + VIR_ALWAYS_INLINE \ + std::enable_if_t, stdx::simd> \ + name (const stdx::simd& x, \ + const vecmath_detail::nondeduced_t>& y) \ + { return std::experimental::parallelism_v2::name(x, y); } \ + \ + template >> \ + VIR_ALWAYS_INLINE stdx::simd \ + name (U&& x, const stdx::simd& y) \ + { \ + /* Qualified, so that argument-dependent lookup cannot add the \ + * underlying overload, which deduces both parameters and would \ + * therefore win partial ordering. */ \ + return vir::vecmath::name(stdx::simd(static_cast(x)), y); \ + } \ + } + +/* Functions the vector math library on this system does not have + * + * They still get a vir::vecmath overload, forwarding to the underlying + * implementation. Whether a function is routed is this header's business, not + * the caller's: code calling vir::vecmath::sinh should compile against any + * glibc and simply be faster on the ones that can. + */ +#define VIR_VECMATH_FN_1_FORWARD(name) \ + namespace vir::vecmath { \ + template \ + VIR_ALWAYS_INLINE \ + std::enable_if_t, stdx::simd> \ + name (const stdx::simd& x) \ + { return std::experimental::parallelism_v2::name(x); } \ + } + +#define VIR_VECMATH_FN_2_FORWARD(name) \ + namespace vir::vecmath { \ + template \ + VIR_ALWAYS_INLINE \ + std::enable_if_t, stdx::simd> \ + name (const stdx::simd& x, \ + const vecmath_detail::nondeduced_t>& y) \ + { return std::experimental::parallelism_v2::name(x, y); } \ + \ + template >> \ + VIR_ALWAYS_INLINE stdx::simd \ + name (U&& x, const stdx::simd& y) \ + { return vir::vecmath::name(stdx::simd(static_cast(x)), y); } \ + } + +// available since glibc 2.22 +VIR_VECMATH_FN_1(sin) +VIR_VECMATH_FN_1(cos) +VIR_VECMATH_FN_1(exp) +VIR_VECMATH_FN_1(log) +VIR_VECMATH_FN_2(pow) + +#ifdef VIR_HAVE_SIMD_VECMATH_EXTENDED +VIR_VECMATH_FN_1(tan) +VIR_VECMATH_FN_1(asin) +VIR_VECMATH_FN_1(acos) +VIR_VECMATH_FN_1(atan) +VIR_VECMATH_FN_1(sinh) +VIR_VECMATH_FN_1(cosh) +VIR_VECMATH_FN_1(tanh) +VIR_VECMATH_FN_1(asinh) +VIR_VECMATH_FN_1(acosh) +VIR_VECMATH_FN_1(atanh) +VIR_VECMATH_FN_1(exp2) +VIR_VECMATH_FN_1(expm1) +VIR_VECMATH_FN_1(log2) +VIR_VECMATH_FN_1(log10) +VIR_VECMATH_FN_1(log1p) +VIR_VECMATH_FN_1(cbrt) +VIR_VECMATH_FN_1(erf) +VIR_VECMATH_FN_1(erfc) +VIR_VECMATH_FN_2(atan2) +#else +VIR_VECMATH_FN_1_FORWARD(tan) +VIR_VECMATH_FN_1_FORWARD(asin) +VIR_VECMATH_FN_1_FORWARD(acos) +VIR_VECMATH_FN_1_FORWARD(atan) +VIR_VECMATH_FN_1_FORWARD(sinh) +VIR_VECMATH_FN_1_FORWARD(cosh) +VIR_VECMATH_FN_1_FORWARD(tanh) +VIR_VECMATH_FN_1_FORWARD(asinh) +VIR_VECMATH_FN_1_FORWARD(acosh) +VIR_VECMATH_FN_1_FORWARD(atanh) +VIR_VECMATH_FN_1_FORWARD(exp2) +VIR_VECMATH_FN_1_FORWARD(expm1) +VIR_VECMATH_FN_1_FORWARD(log2) +VIR_VECMATH_FN_1_FORWARD(log10) +VIR_VECMATH_FN_1_FORWARD(log1p) +VIR_VECMATH_FN_1_FORWARD(cbrt) +VIR_VECMATH_FN_1_FORWARD(erf) +VIR_VECMATH_FN_1_FORWARD(erfc) +VIR_VECMATH_FN_2_FORWARD(atan2) +#endif + +/* Deliberately not routed here: + * + * hypot, because the underlying implementation already evaluates it with SIMD + * instructions, including the fixups libmvec's version would not do, and + * because its overload set (two and three arguments, plus the converting + * forms) is larger than the pair generated above. + */ + +#endif // VIR_HAVE_SIMD_VECMATH +#endif // VIR_SIMD_VECMATH_H_