diff --git a/Common/Math/CrossSIMD.h b/Common/Math/CrossSIMD.h
index 18a7526d47..d66962333a 100644
--- a/Common/Math/CrossSIMD.h
+++ b/Common/Math/CrossSIMD.h
@@ -159,6 +159,9 @@ struct Vec4S32 {
 		// in the other register.
 		return Vec4S32{ _mm_madd_epi16(v, _mm_and_si128(other.v, _mm_set1_epi32(0x0000FFFF))) };
 	}
+	// Each 16-bit half times the other's, keeping the low 16 bits of each product: two lanes of 16 bits
+	// in each 32-bit one.
+	Vec4S32 MulHalves(Vec4S32 other) const { return Vec4S32{ _mm_mullo_epi16(v, other.v) }; }
 
 	Vec4S32 SignExtend16() const { return Vec4S32{ _mm_srai_epi32(_mm_slli_epi32(v, 16), 16) }; }
 	// NOTE: These can be done in sequence, but when done, you must FixupAfterMinMax to get valid output.
@@ -647,6 +650,9 @@ struct Vec4S32 {
 
 	// Warning: Unlike on x86, this is a full 32-bit multiplication.
 	Vec4S32 Mul16(Vec4S32 other) const { return Vec4S32{ vmulq_s32(v, other.v) }; }
+	// Each 16-bit half times the other's, keeping the low 16 bits of each product: two lanes of 16 bits
+	// in each 32-bit one.
+	Vec4S32 MulHalves(Vec4S32 other) const { return Vec4S32{ vreinterpretq_s32_s16(vmulq_s16(vreinterpretq_s16_s32(v), vreinterpretq_s16_s32(other.v))) }; }
 
 	Vec4S32 SignExtend16() const { return Vec4S32{ vshrq_n_s32(vshlq_n_s32(v, 16), 16) }; }
 	// NOTE: These can be done in sequence, but when done, you must FixupAfterMinMax to get valid output (on SSE2 at least).
@@ -1214,6 +1220,9 @@ struct Vec4S32 {
 
 	// Warning: Unlike on x86, this is a full 32-bit multiplication.
 	Vec4S32 Mul16(Vec4S32 other) const { return Vec4S32{ __lsx_vmul_w(v, other.v) }; }
+	// Each 16-bit half times the other's, keeping the low 16 bits of each product: two lanes of 16 bits
+	// in each 32-bit one.
+	Vec4S32 MulHalves(Vec4S32 other) const { return Vec4S32{ __lsx_vmul_h(v, other.v) }; }
 
 	Vec4S32 SignExtend16() const { return Vec4S32{ __lsx_vsrai_w(__lsx_vslli_w(v, 16), 16) }; }
 	// NOTE: These can be done in sequence, but when done, you must FixupAfterMinMax to get valid output (on SSE2 at least).
@@ -1696,6 +1705,17 @@ struct Vec4S32 {
 
 	// Warning: Unlike on x86 SSE2, this is a full 32-bit multiplication.
 	Vec4S32 Mul16(Vec4S32 other) const { return Vec4S32{ { v[0] * other.v[0], v[1] * other.v[1], v[2] * other.v[2], v[3] * other.v[3] } }; }
+	// Each 16-bit half times the other's, keeping the low 16 bits of each product: two lanes of 16 bits
+	// in each 32-bit one.
+	Vec4S32 MulHalves(Vec4S32 other) const {
+		Vec4S32 r;
+		for (int i = 0; i < 4; ++i) {
+			const uint32_t lo = ((uint32_t)v[i] * (uint32_t)other.v[i]) & 0xFFFF;
+			const uint32_t hi = ((uint32_t)v[i] >> 16) * ((uint32_t)other.v[i] >> 16);
+			r.v[i] = (int32_t)((hi << 16) | lo);
+		}
+		return r;
+	}
 
 	Vec4S32 SignExtend16() const {
 		Vec4S32 tmp;
diff --git a/GPU/Software/DrawPixel.cpp b/GPU/Software/DrawPixel.cpp
index b618c49913..ac873e4ca4 100644
--- a/GPU/Software/DrawPixel.cpp
+++ b/GPU/Software/DrawPixel.cpp
@@ -938,9 +938,9 @@ static inline void DrawSpanVector(int x, int y, const int *maskIn, const int *zI
 			// Like the BLEND texfunc, always rounding up.
 			const Vec4S32 f = Vec4S32::Load(fogIn), invF = c255 - f;
 			const uint32_t fc = pixelID.cached.fogColor;
-			r = (r.Mul(f) + Vec4S32::Splat(fc & 0xFF).Mul(invF) + c255).Shr<8>();
-			g = (g.Mul(f) + Vec4S32::Splat((fc >> 8) & 0xFF).Mul(invF) + c255).Shr<8>();
-			b = (b.Mul(f) + Vec4S32::Splat((fc >> 16) & 0xFF).Mul(invF) + c255).Shr<8>();
+			r = (r.Mul16(f) + Vec4S32::Splat(fc & 0xFF).Mul16(invF) + c255).Shr<8>();
+			g = (g.Mul16(f) + Vec4S32::Splat((fc >> 8) & 0xFF).Mul16(invF) + c255).Shr<8>();
+			b = (b.Mul16(f) + Vec4S32::Splat((fc >> 16) & 0xFF).Mul16(invF) + c255).Shr<8>();
 		}
 
 		if (pixelID.colorTest) {
@@ -1021,7 +1021,7 @@ static inline void DrawSpanVector(int x, int y, const int *maskIn, const int *zI
 	if constexpr (blend == SpanBlend::SRC_ALPHA) {
 		const Vec4S32 one = Vec4S32::Splat(1);
 		auto term = [&](Vec4S32 v, Vec4S32 f) {
-			return (v.Shl<1>() + one).Mul(f.Shl<1>() + one).Shr<10>();
+			return (v.Shl<1>() + one).Mul16(f.Shl<1>() + one).Shr<10>();
 		};
 		const Vec4S32 invA = c255 - a;
 		r = term(r, a) + term(dst[0], invA);
@@ -1057,7 +1057,7 @@ static inline void DrawSpanVector(int x, int y, const int *maskIn, const int *zI
 		Vec4S32 out[3];
 		for (int c = 0; c < 3; ++c) {
 			auto term = [&](Vec4S32 v, Vec4S32 f) {
-				return (v.Shl<1>() + one).Mul(f.Shl<1>() + one).Shr<10>();
+				return (v.Shl<1>() + one).Mul16(f.Shl<1>() + one).Shr<10>();
 			};
 			switch (pixelID.AlphaBlendEq()) {
 			case GE_BLENDMODE_MUL_AND_ADD: out[c] = term(src[c], sf[c]) + term(dst[c], df[c]); break;
diff --git a/GPU/Software/Sampler.cpp b/GPU/Software/Sampler.cpp
index 68a3b6f9f9..04d4afc9f1 100644
--- a/GPU/Software/Sampler.cpp
+++ b/GPU/Software/Sampler.cpp
@@ -710,12 +710,17 @@ static inline void TexelPairs4(Vec4F32 st, int size, bool clamp, Vec4S32 &c0, Ve
 	}
 }
 
-// RowOffsetT and ColumnOffsetT for four texels at once.
+// v times k, both from 0 to 32767 for Mul16 (a 32-bit multiply is slow on SSE2).
+static inline Vec4S32 MulRow(Vec4S32 v, int k) {
+	return k < 32768 ? v.Mul16(Vec4S32::Splat(k)) : v.Mul(Vec4S32::Splat(k));
+}
+
+// RowOffsetT and ColumnOffsetT for four texels at once. The rows are below 512.
 template <uint32_t bits, bool swizzled>
 static inline Vec4S32 RowOffsets4(int bufw, Vec4S32 v) {
 	if (!swizzled)
-		return v.Mul(Vec4S32::Splat(bufw * bits >> 3));
-	return v.Shr<3>().Mul(Vec4S32::Splat((bufw * bits / 32) * 32)) + (v & Vec4S32::Splat(7)).Shl<4>();
+		return MulRow(v, bufw * bits >> 3);
+	return MulRow(v.Shr<3>(), bufw * bits / 32).Shl<5>() + (v & Vec4S32::Splat(7)).Shl<4>();
 }
 
 template <uint32_t bits, bool swizzled>
@@ -791,9 +796,14 @@ static inline void SampleLinearLevel4LanesT(Vec4F32 s, Vec4F32 t, const u8 *tptr
 	const Vec4S32 tl = Vec4S32::LoadAligned((const int *)texels[0]), tr = Vec4S32::LoadAligned((const int *)texels[1]);
 	const Vec4S32 bl = Vec4S32::LoadAligned((const int *)texels[2]), br = Vec4S32::LoadAligned((const int *)texels[3]);
 	const Vec4S32 halves = Vec4S32::Splat(0x00FF00FF);
-	const Vec4S32 invU = Vec4S32::Splat(16) - fracU, invV = Vec4S32::Splat(16) - fracV;
+	// The fractions in both halves, for MulHalves: no product exceeds 16 bits either.
+	const Vec4S32 sixteen = Vec4S32::Splat(0x00100010);
+	const Vec4S32 fracU2 = fracU | fracU.Shl<16>(), fracV2 = fracV | fracV.Shl<16>();
+	const Vec4S32 invU = sixteen - fracU2, invV = sixteen - fracV2;
+	fracU = fracU2;
+	fracV = fracV2;
 	auto lerp = [&](Vec4S32 a, Vec4S32 b, Vec4S32 inv, Vec4S32 f) {
-		return (a.Mul(inv) + b.Mul(f)).Shr<4>() & halves;
+		return (a.MulHalves(inv) + b.MulHalves(f)).Shr<4>() & halves;
 	};
 	Vec4S32 pairs[2];
 	for (int h = 0; h < 2; ++h) {
@@ -819,11 +829,11 @@ static inline void TextureFunction4(const Vec4S32 prim[4], const Vec4S32 tex[4],
 	const bool doubling = samplerID.useColorDoubling;
 	const Vec4S32 one = Vec4S32::Splat(1), c255 = Vec4S32::Splat(255);
 	// The alpha blended the common way: (prim + 1) * tex / 256.
-	const Vec4S32 modulatedA = rgba ? (prim[3] + one).Mul(tex[3]).Shr<8>() : prim[3];
+	const Vec4S32 modulatedA = rgba ? (prim[3] + one).Mul16(tex[3]).Shr<8>() : prim[3];
 	switch (samplerID.TexFunc()) {
 	case GE_TEXFUNC_MODULATE:
 		for (int c = 0; c < 3; ++c)
-			out[c] = (prim[c] + one).Mul(doubling ? tex[c].Shl<1>() : tex[c]).Shr<8>();
+			out[c] = (prim[c] + one).Mul16(doubling ? tex[c].Shl<1>() : tex[c]).Shr<8>();
 		out[3] = modulatedA;
 		break;
 
@@ -832,7 +842,7 @@ static inline void TextureFunction4(const Vec4S32 prim[4], const Vec4S32 tex[4],
 			// Both colors are boosted here, making the alpha have more weight.
 			const Vec4S32 t = tex[3], invt = c255 - tex[3];
 			for (int c = 0; c < 3; ++c) {
-				const Vec4S32 sum = (prim[c] + one).Mul(invt) + (tex[c] + one).Mul(t);
+				const Vec4S32 sum = (prim[c] + one).Mul16(invt) + (tex[c] + one).Mul16(t);
 				out[c] = doubling ? sum.Shr<7>() : sum.Shr<8>();
 			}
 		} else {
@@ -847,7 +857,7 @@ static inline void TextureFunction4(const Vec4S32 prim[4], const Vec4S32 tex[4],
 		// Unlike the others (and even alpha), this one always rounds up.
 		const uint32_t env = samplerID.cached.texBlendColor;
 		for (int c = 0; c < 3; ++c) {
-			const Vec4S32 sum = (c255 - tex[c]).Mul(prim[c]) + tex[c].Mul(Vec4S32::Splat((env >> (8 * c)) & 0xFF)) + c255;
+			const Vec4S32 sum = (c255 - tex[c]).Mul16(prim[c]) + tex[c].Mul16(Vec4S32::Splat((env >> (8 * c)) & 0xFF)) + c255;
 			out[c] = doubling ? sum.Shr<7>() : sum.Shr<8>();
 		}
 		out[3] = modulatedA;
diff --git a/unittest/TestCrossSIMD.cpp b/unittest/TestCrossSIMD.cpp
index 56dd66e382..aa648abfc4 100644
--- a/unittest/TestCrossSIMD.cpp
+++ b/unittest/TestCrossSIMD.cpp
@@ -445,6 +445,19 @@ static bool TestVec4F32Loads() {
 	EXPECT_EQ_FLOAT(result[2], 1e20f);
 	EXPECT_EQ_FLOAT(result[3], 1.0f);
 
+	// MulHalves: two 16-bit products per lane, each kept to 16 bits.
+	{
+		const int halvesA[4] = { 0x00FF00FF, 0x01230045, (int)0xFFFF0002, 0x7FFF8000 };
+		const int halvesB[4] = { 0x00100010, 0x00030007, 0x00020003, 0x00020002 };
+		int halvesOut[4];
+		Vec4S32::Load(halvesA).MulHalves(Vec4S32::Load(halvesB)).Store(halvesOut);
+		for (int i = 0; i < 4; i++) {
+			const uint32_t lo = ((uint32_t)halvesA[i] * (uint32_t)halvesB[i]) & 0xFFFF;
+			const uint32_t hi = (((uint32_t)halvesA[i] >> 16) * ((uint32_t)halvesB[i] >> 16)) & 0xFFFF;
+			EXPECT_EQ_HEX((uint32_t)halvesOut[i], (hi << 16) | lo);
+		}
+	}
+
 	// LoadS64Low: the low halves.
 	int iresult[4];
 	const int64_t wide[4] = { 0x123456789LL, -2, (int64_t)0x7FFFFFFF << 1, -0x100000000LL + 7 };