From fe2c9a3231c8d338be31cd0dd34ea406364f1f9b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Henrik=20Rydg=C3=A5rd?= Date: Wed, 14 Oct 2026 23:11:12 +0200 Subject: [PATCH] VFPU: Add scalar and LSX versions of vfpu_dot Hosts without NEON or SSE2 ran vfpu_dot_reference() for every dot product. vfpu_dot_scalar() is the same computation without its branches, finished like the SIMD versions, and is what RISC-V gets. LoongArch builds with LSX (GCC's default la64v1.0 has it) get vfpu_dot_lsx(), the SSE2 version in LSX. VFPUDot checks vfpu_dot_scalar() on every host, and its inputs come from ForEachVFPUDotTestCase(), for the JITs' routines to use too. Co-Authored-By: Claude Opus 5.5 (1M context) --- Core/MIPS/MIPSVFPUUtils.cpp | 78 ++++++++++++++++++++++++++++++++++--- Core/MIPS/MIPSVFPUUtils.h | 2 + unittest/UnitTest.cpp | 54 ++++++++++++++++++------- unittest/UnitTest.h | 5 +++ 4 files changed, 120 insertions(+), 19 deletions(-) diff --git a/Core/MIPS/MIPSVFPUUtils.cpp b/Core/MIPS/MIPSVFPUUtils.cpp index 6c55d2e518..58c1be3229 100644 --- a/Core/MIPS/MIPSVFPUUtils.cpp +++ b/Core/MIPS/MIPSVFPUUtils.cpp @@ -867,14 +867,12 @@ float vfpu_dot_reference(const float a[4], const float b[4]) { return ret; } -#if PPSSPP_ARCH(ARM64_NEON) || PPSSPP_ARCH(SSE2) - // Inf and NaN inputs are rare, and the reference handles them. static NO_INLINE float vfpu_dot_special(const float a[4], const float b[4]) { return vfpu_dot_reference(a, b); } -// The end of the SIMD versions: drops the extra bits off the exact sum by truncation, then rounds to +// The end of the fast versions: drops the extra bits off the exact sum by truncation, then rounds to // 24 bits, nearest even. With the top bit moved to bit 62, adding 2^38 - 1 plus the lowest kept bit // and shifting out 39 bits does that for every magnitude, and adding the significand, implicit bit // and all, onto the exponent field lets a carry from the rounding bump the exponent by itself. It's @@ -894,7 +892,35 @@ static inline float vfpu_dot_finish(int32_t val, uint32_t ehi) { return ret; } -#endif +// The reference without its branches, for hosts with no SIMD version. Bit-exact with it. +float vfpu_dot_scalar(const float a[4], const float b[4]) { + uint32_t x[4], y[4], e[4]; + memcpy(x, a, sizeof(x)); + memcpy(y, b, sizeof(y)); + + // Exponent sums, zero for zero and subnormal inputs, and whether any input is inf or NaN. + uint32_t ehi = 0, special = 0; + for (int i = 0; i < 4; i++) { + const uint32_t ex = (x[i] >> 23) & 0xFF, ey = (y[i] >> 23) & 0xFF; + special |= (ex + 1) | (ey + 1); + e[i] = ex != 0 && ey != 0 ? ex + ey : 0; + ehi = e[i] > ehi ? e[i] : ehi; + } + if (special & 0x100) + return vfpu_dot_special(a, b); + + // 24x24-bit products kept to 2 extra bits with round-to-odd, aligned to the largest exponent + // by truncation (p is below 2^27, so 31 does for any larger shift), and summed exactly. + int32_t val = 0; + for (int i = 0; i < 4; i++) { + const uint64_t v = (uint64_t)((x[i] & 0x007FFFFF) | 0x00800000) * ((y[i] & 0x007FFFFF) | 0x00800000); + uint32_t p = (uint32_t)(v >> 21) | ((v & ((1 << 21) - 1)) != 0); + const uint32_t d = ehi - e[i]; + p = e[i] != 0 ? p >> (d < 31 ? d : 31) : 0; + val += (int32_t)((x[i] ^ y[i]) >> 31 ? 0u - p : p); + } + return vfpu_dot_finish(val, ehi); +} #if PPSSPP_ARCH(ARM64_NEON) @@ -983,6 +1009,46 @@ static float vfpu_dot_sse2(const float a[4], const float b[4]) { return vfpu_dot_finish(_mm_cvtsi128_si32(v), ehi); } +#elif PPSSPP_ARCH(LOONGARCH64_LSX) + +// The reference, four lanes at a time, the same way as vfpu_dot_sse2(). Bit-exact with it. +static float vfpu_dot_lsx(const float a[4], const float b[4]) { + const __m128i x = __lsx_vld(a, 0); + const __m128i y = __lsx_vld(b, 0); + const __m128i expBits = __lsx_vreplgr2vr_w(0x7F800000); + const __m128i xe = __lsx_vand_v(x, expBits), ye = __lsx_vand_v(y, expBits); + + // Zero and subnormal inputs make the product zero, with the lowest exponent. Inf and NaN get + // 0x7FF, above any real exponent sum, so one maximum finds the alignment and the specials. + const __m128i flush = __lsx_vor_v(__lsx_vseqi_w(xe, 0), __lsx_vseqi_w(ye, 0)); + const __m128i special = __lsx_vor_v(__lsx_vseq_w(xe, expBits), __lsx_vseq_w(ye, expBits)); + const __m128i e = __lsx_vor_v(__lsx_vandn_v(flush, __lsx_vsrli_w(__lsx_vadd_w(xe, ye), 23)), __lsx_vsrli_w(special, 21)); + __m128i emax = __lsx_vmax_wu(e, __lsx_vshuf4i_w(e, 0x4E)); + emax = __lsx_vmax_wu(emax, __lsx_vshuf4i_w(emax, 0xB1)); + const uint32_t ehi = (uint32_t)__lsx_vpickve2gr_w(emax, 0); + if (ehi >= 0x7FF) + return vfpu_dot_special(a, b); + + // 24x24-bit products, even and odd lanes separately, kept to 2 extra bits with round-to-odd. + const __m128i mantMask = __lsx_vreplgr2vr_w(0x007FFFFF), implicitBit = __lsx_vreplgr2vr_w(0x00800000); + const __m128i mx = __lsx_vor_v(__lsx_vand_v(x, mantMask), implicitBit); + const __m128i my = __lsx_vor_v(__lsx_vand_v(y, mantMask), implicitBit); + const __m128i pEven = __lsx_vmulwev_d_wu(mx, my); + const __m128i pOdd = __lsx_vmulwod_d_wu(mx, my); + const __m128i kept = __lsx_vpackev_w(__lsx_vsrli_d(pOdd, 21), __lsx_vsrli_d(pEven, 21)); + const __m128i dropped = __lsx_vand_v(__lsx_vpackev_w(pOdd, pEven), __lsx_vreplgr2vr_w((1 << 21) - 1)); + const __m128i p = __lsx_vandn_v(flush, __lsx_vor_v(kept, __lsx_vmini_wu(dropped, 1))); + + // Align to the largest exponent by truncation (p is below 2^27, so 31 does for any larger + // shift), apply the signs, and sum exactly. + const __m128i d = __lsx_vmini_wu(__lsx_vsub_w(emax, e), 31); + const __m128i sign = __lsx_vsrai_w(__lsx_vxor_v(x, y), 31); + __m128i v = __lsx_vsub_w(__lsx_vxor_v(__lsx_vsrl_w(p, d), sign), sign); + v = __lsx_vadd_w(v, __lsx_vshuf4i_w(v, 0x4E)); + v = __lsx_vadd_w(v, __lsx_vshuf4i_w(v, 0xB1)); + return vfpu_dot_finish(__lsx_vpickve2gr_w(v, 0), ehi); +} + #endif float vfpu_dot(const float a[4], const float b[4]) { @@ -990,8 +1056,10 @@ float vfpu_dot(const float a[4], const float b[4]) { return vfpu_dot_neon(a, b); #elif PPSSPP_ARCH(SSE2) return vfpu_dot_sse2(a, b); +#elif PPSSPP_ARCH(LOONGARCH64_LSX) + return vfpu_dot_lsx(a, b); #else - return vfpu_dot_reference(a, b); + return vfpu_dot_scalar(a, b); #endif } diff --git a/Core/MIPS/MIPSVFPUUtils.h b/Core/MIPS/MIPSVFPUUtils.h index e3eb9572b0..ac5f603344 100644 --- a/Core/MIPS/MIPSVFPUUtils.h +++ b/Core/MIPS/MIPSVFPUUtils.h @@ -74,6 +74,8 @@ inline float vfpu_clamp(float v, float min, float max) { float vfpu_dot(const float a[4], const float b[4]); // The portable version vfpu_dot is checked against. float vfpu_dot_reference(const float a[4], const float b[4]); +// vfpu_dot on hosts with no SIMD version. +float vfpu_dot_scalar(const float a[4], const float b[4]); // Four vfpu_dot()s against v, over the first terms lanes: out[k] uses block[4 * i + k] for lane i, // or block[4 * k + i] if transposed (see IROp::Vec4VfpuMultiDot). void vfpu_multi_dot(float out[4], const float block[16], const float v[4], int terms, bool transposed); diff --git a/unittest/UnitTest.cpp b/unittest/UnitTest.cpp index 1bc4a63fb4..16af0e29b5 100644 --- a/unittest/UnitTest.cpp +++ b/unittest/UnitTest.cpp @@ -1823,10 +1823,10 @@ bool TestFastVec() { return true; } -// vfpu_dot's SIMD versions against the reference, on inputs chosen to make trouble: close -// exponents, cancelling products, ties, zeroes and subnormals, the overflow and underflow edges, -// inf and NaN, and sums whose rounding carries into the next power of two. -bool TestVFPUDot() { +// Calls check on inputs chosen to make trouble for a VFPU dot product: close exponents, cancelling +// products, ties, zeroes and subnormals, the overflow and underflow edges, inf and NaN, and sums whose +// rounding carries into the next power of two. Stops at the first false. +bool ForEachVFPUDotTestCase(const std::function &check) { uint64_t state = 0x9E3779B97F4A7C15ULL; auto rnd = [&]() { state ^= state << 13; @@ -1844,16 +1844,6 @@ bool TestVFPUDot() { memcpy(&bits, &f, sizeof(bits)); return bits; }; - auto check = [&](const float a[4], const float b[4]) { - const uint32_t expected = toBits(vfpu_dot_reference(a, b)); - const uint32_t actual = toBits(vfpu_dot(a, b)); - if (expected != actual) { - printf("vfpu_dot(%08x %08x %08x %08x, %08x %08x %08x %08x) = %08x, expected %08x\n", - toBits(a[0]), toBits(a[1]), toBits(a[2]), toBits(a[3]), toBits(b[0]), toBits(b[1]), toBits(b[2]), toBits(b[3]), actual, expected); - return false; - } - return true; - }; for (int n = 0; n < 4000000; n++) { float a[4], b[4]; const int mode = (int)(rnd() & 15); @@ -1891,6 +1881,20 @@ bool TestVFPUDot() { return false; } + // Cases that went wrong somewhere. + static const uint32_t fixed[][8] = { + { 0x8055f97b, 0xcaa86851, 0, 0, 0x44282000, 0x3f800000, 0, 0 }, + }; + for (const auto &f : fixed) { + float a[4], b[4]; + for (int i = 0; i < 4; i++) { + a[i] = fromBits(f[i]); + b[i] = fromBits(f[i + 4]); + } + if (!check(a, b)) + return false; + } + // Sums just below and above a power of two, whose rounding carries into the next exponent: // 1.0 from just under it, and inf at the top of the range. for (int e = 1; e <= 254; e++) { @@ -1915,6 +1919,28 @@ bool TestVFPUDot() { return true; } + +// vfpu_dot's fast versions against the reference. +bool TestVFPUDot() { + auto toBits = [](float f) { + uint32_t bits; + memcpy(&bits, &f, sizeof(bits)); + return bits; + }; + auto check = [&](const float a[4], const float b[4]) { + const uint32_t expected = toBits(vfpu_dot_reference(a, b)); + const uint32_t actual = toBits(vfpu_dot(a, b)); + const uint32_t scalar = toBits(vfpu_dot_scalar(a, b)); + if (expected != actual || expected != scalar) { + printf("vfpu_dot(%08x %08x %08x %08x, %08x %08x %08x %08x) = %08x (scalar %08x), expected %08x\n", + toBits(a[0]), toBits(a[1]), toBits(a[2]), toBits(a[3]), toBits(b[0]), toBits(b[1]), toBits(b[2]), toBits(b[3]), actual, scalar, expected); + return false; + } + return true; + }; + return ForEachVFPUDotTestCase(check); +} + bool TestVFPUSinCos() { float sine, cosine; vfpu_sincos(0.0f, sine, cosine); diff --git a/unittest/UnitTest.h b/unittest/UnitTest.h index 1d38e94a2d..963e43609d 100644 --- a/unittest/UnitTest.h +++ b/unittest/UnitTest.h @@ -4,6 +4,7 @@ #include #include #include +#include #include "Common/TimeUtil.h" @@ -48,3 +49,7 @@ inline bool rel_equal(float a, float b, float precision) { #define RET(a) if (!(a)) { return false; } extern bool g_testLog; + +// Calls check on a few million inputs chosen to make trouble for a VFPU dot product (see +// UnitTest.cpp), and stops at the first false. +bool ForEachVFPUDotTestCase(const std::function &check);