Add tests on transcendental values

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<float>::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

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<float>::infinity();
+	const float nan = std::numeric_limits<float>::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<float>::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<int W>
+void rsqrtRcpRangeTest()
+{
+	using weaver::all;
+	using weaver::equal;
+
+	using floatN = weaver::VectorScalar<float, W>;
+
+	// 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;