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) <noreply@anthropic.com>

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
 }
 

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);

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<bool(const float a[4], const float b[4])> &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);

unittest/UnitTest.h

@@ -4,6 +4,7 @@
 #include <cstdio>
 #include <cstring>
 #include <algorithm>
+#include <functional>
 
 #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<bool(const float a[4], const float b[4])> &check);