From 8cfa3982ea937d0b484a6175570bb0734e71404a Mon Sep 17 00:00:00 2001 From: Basil Fierz Date: Sat, 12 Sep 2026 00:14:08 +0200 Subject: [PATCH] Add tests on transcendental values --- src/intrinsics_neon.cpp | 8 +- src/tests/simd.cpp | 183 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 190 insertions(+), 1 deletion(-) diff --git a/src/intrinsics_neon.cpp b/src/intrinsics_neon.cpp index 3b0c9d9..efcb8a3 100644 --- a/src/intrinsics_neon.cpp +++ b/src/intrinsics_neon.cpp @@ -123,7 +123,13 @@ namespace weaver { float32x4_t weaver_vsqrtq_f32(float32x4_t x) { const float4 v{ x }; - return select(v == 0.0f, 0.0f, v * rsqrt(v)).get(0); + const float4 posInf{ std::numeric_limits::infinity() }; + + // v * rsqrt(v) breaks down for +Inf: the rsqrt estimate underflows to + // 0 for such a large input, so the product becomes Inf * 0 = NaN + // instead of the correct +Inf. Zero is special-cased for the + // analogous reason (0 * Inf = NaN). + return select(v == posInf, posInf, select(v == 0.0f, float4(0.0f), v * rsqrt(v))).get(0); } } #endif diff --git a/src/tests/simd.cpp b/src/tests/simd.cpp index 7c6ca10..e15df43 100644 --- a/src/tests/simd.cpp +++ b/src/tests/simd.cpp @@ -217,6 +217,189 @@ TEST_CASE("Log FloatN", "[float][operation]") CHECK(all(equal(ref3, res3, float16(1e-5f)))); } +TEST_CASE("Sqrt FloatN special values", "[float][operation]") +{ + using weaver::all; + using weaver::equal; + using weaver::float16; + using weaver::float4; + using weaver::float8; + using weaver::isinf; + + const float inf = std::numeric_limits::infinity(); + const float nan = std::numeric_limits::quiet_NaN(); + + SECTION("Negative input yields NaN") + { + float4 res1 = sqrt(float4(-1.0f)); + float8 res2 = sqrt(float8(-1.0f)); + float16 res3 = sqrt(float16(-1.0f)); + + for (int i = 0; i < 4; i++) CHECK(std::isnan(res1[i])); + for (int i = 0; i < 8; i++) CHECK(std::isnan(res2[i])); + for (int i = 0; i < 16; i++) CHECK(std::isnan(res3[i])); + } + + SECTION("NaN input propagates") + { + float4 res1 = sqrt(float4(nan)); + float8 res2 = sqrt(float8(nan)); + float16 res3 = sqrt(float16(nan)); + + for (int i = 0; i < 4; i++) CHECK(std::isnan(res1[i])); + for (int i = 0; i < 8; i++) CHECK(std::isnan(res2[i])); + for (int i = 0; i < 16; i++) CHECK(std::isnan(res3[i])); + } + + SECTION("Infinity input yields infinity") + { + CHECK(all(isinf(sqrt(float4(inf))))); + CHECK(all(isinf(sqrt(float8(inf))))); + CHECK(all(isinf(sqrt(float16(inf))))); + } + + SECTION("Zero input yields zero") + { + CHECK(all(equal(sqrt(float4(0.0f)), float4(0.0f)))); + CHECK(all(equal(sqrt(float8(0.0f)), float8(0.0f)))); + CHECK(all(equal(sqrt(float16(0.0f)), float16(0.0f)))); + } + + SECTION("Denormal input matches std::sqrt") + { + const float denorm = std::numeric_limits::denorm_min(); + const float ref = std::sqrt(denorm); + + CHECK(all(equal(sqrt(float4(denorm)), float4(ref), float4(1e-5f)))); + CHECK(all(equal(sqrt(float8(denorm)), float8(ref), float8(1e-5f)))); + CHECK(all(equal(sqrt(float16(denorm)), float16(ref), float16(1e-5f)))); + } +} + +template +void rsqrtRcpRangeTest() +{ + using weaver::all; + using weaver::equal; + + using floatN = weaver::VectorScalar; + + // Values spanning many orders of magnitude, but staying well clear of the + // range near FLT_MAX where the Newton-Raphson refinement's intermediate + // nr*nr term can overflow and the approximation breaks down. + const float values[] = { 1e-10f, 1e-6f, 1e-3f, 1.0f, 7.3f, 1e3f, 1e6f, 1e10f, 1e15f }; + + for (float v : values) + { + const floatN vec(v); + const float refRsqrt = 1.0f / std::sqrt(v); + const float refRcp = 1.0f / v; + + CAPTURE(v); + CHECK(all(equal(rsqrt(vec), floatN(refRsqrt), floatN(1e-4f)))); + CHECK(all(equal(rcp(vec), floatN(refRcp), floatN(1e-4f)))); + } +} + +TEST_CASE("Rsqrt/Rcp FloatN across magnitude range", "[float][operation]") +{ + rsqrtRcpRangeTest<4>(); + rsqrtRcpRangeTest<8>(); + rsqrtRcpRangeTest<16>(); +} + +TEST_CASE("Rcp FloatN zero input", "[float][operation]") +{ + using weaver::all; + using weaver::float16; + using weaver::float4; + using weaver::float8; + using weaver::isinf; + + // rcp(0) is explicitly special-cased to +Inf rather than propagating the + // NaN that 0 * Inf would otherwise produce in the Newton-Raphson refinement. + CHECK(all(isinf(rcp(float4(0.0f))))); + CHECK(all(isinf(rcp(float8(0.0f))))); + CHECK(all(isinf(rcp(float16(0.0f))))); +} + +TEST_CASE("Log FloatN boundary values", "[float][operation]") +{ + using weaver::all; + using weaver::equal; + using weaver::float16; + using weaver::float4; + using weaver::float8; + + SECTION("Non-positive input is non-finite") + { + // The vectorized (SSE/AVX/AVX512/NEON) log_ps is documented to return + // NaN for x <= 0; the ref backend forwards to std::log, which returns + // -Inf for 0 and NaN for negative input. Only assert what holds for + // both: negative input is NaN, and zero is at least non-finite. + float4 res1 = log(float4(0.0f)); + float4 res2 = log(float4(-1.0f)); + + for (int i = 0; i < 4; i++) + { + CHECK_FALSE(std::isfinite(res1[i])); + CHECK(std::isnan(res2[i])); + } + } + + SECTION("Values straddling the internal sqrt(1/2) branch threshold") + { + // The polynomial approximation switches behavior around SQRTHF + // (~0.70710678), so values just below/above it are the most likely + // place for the approximation to diverge from std::log. + const float values[] = { 0.7071065f, 0.70710678f, 0.7071072f }; + + for (float v : values) + { + const float ref = std::log(v); + + CAPTURE(v); + CHECK(all(equal(log(float4(v)), float4(ref), float4(1e-5f)))); + CHECK(all(equal(log(float8(v)), float8(ref), float8(1e-5f)))); + CHECK(all(equal(log(float16(v)), float16(ref), float16(1e-5f)))); + } + } + + SECTION("Values straddling one") + { + const float values[] = { 1.0f - 1e-6f, 1.0f, 1.0f + 1e-6f }; + + for (float v : values) + { + const float ref = std::log(v); + + CAPTURE(v); + CHECK(all(equal(log(float4(v)), float4(ref), float4(1e-5f)))); + } + } + + SECTION("Large magnitude input stays finite and close to std::log") + { + const float v = 1e30f; + const float ref = std::log(v); + + CHECK(all(equal(log(float4(v)), float4(ref), float4(1e-4f)))); + } +} + +// Known limitation (not covered by an assertion here because the actual +// outcome turns out to be compiler/codegen-dependent rather than a portable +// contract): the vectorized (SSE/AVX/AVX512/NEON) log/exp implementations +// are ported from *_mathfun.h, which only special-cases x <= 0 for log via a +// plain comparison; that comparison is false for NaN, so on at least some +// compilers a NaN lane is treated as an ordinary (garbage) finite value +// instead of producing NaN. exp has the same gap (its range clamp uses +// min/max, which can silently discard NaN), so pow(x, y) = exp(log(x) * y) +// can likewise fail to propagate NaN, and return a finite garbage value +// instead of 0 for pow(0, y) because log(0) is forced to NaN internally. +// The ref backend forwards straight to std::log/std::exp and does not have +// this problem. + TEST_CASE("Abs FloatN", "[float][operation]") { using weaver::all;