Split MulHigh extensions into CPU and GPU variants

PerfectNumbers.Core/UInt128Extensions.cs

@@ -760,7 +760,7 @@ public static class UInt128Extensions
     // expression tree around the shift. Mirroring the MulHigh layout lets RyuJIT
     // keep the intermediate sum in registers instead of reloading the partial
     // product from the stack.
-    ulong high = aLow.MulHigh(bLow);
+    ulong high = aLow.MulHighCpu(bLow);
     high += aLow * bHigh;
 
     return ((UInt128)high << 64) | low;

PerfectNumbers.Core/ULongExtensions.Gpu.cs

@@ -11,6 +11,7 @@ namespace PerfectNumbers.Core;
 public static partial class ULongExtensions
 {
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	[Obsolete("Use ULongExtensions.MulHighCpu or ULongExtensions.MulHighGpu.")]
 	public static ulong MulHighGpuCompatible(this ulong x, ulong y)
 	{
 		GpuUInt128 product = new(x);

PerfectNumbers.Core/ULongExtensions.cs

@@ -430,14 +430,19 @@ public static partial class ULongExtensions
 	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b);
+	public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHighCpu(b) << 64) | (UInt128)(a * b);
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public static ulong MulHigh(this ulong x, ulong y)
+	public static ulong MulHighCpu(this ulong x, ulong y)
 	{
-		// TODO: Investigate replacing this manual decomposition with the UInt128-based implementation
-		// for CPU callers; the latest benchmarks show the intrinsic path is an order of magnitude
-		// faster, while GPU code can keep using GpuUInt128.MulHigh.
+		return (ulong)(((UInt128)x * y) >> 64);
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong MulHighGpu(this ulong x, ulong y)
+	{
+		// The CPU helpers rely on UInt128 via MulHighCpu; retain this manual decomposition for GPU-style
+		// arithmetic so callers can avoid the intrinsic when staging work for accelerators.
 		ulong xLow = (uint)x;
 		ulong xHigh = x >> 32;
 		ulong yLow = (uint)y;
@@ -554,7 +559,7 @@ public static partial class ULongExtensions
 		ulong m = unchecked(tLow * nPrime);
 		ulong mTimesModulusLow = unchecked(m * modulus);
 
-		ulong result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL));
+		ulong result = unchecked(a.MulHighCpu(b) + m.MulHighCpu(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL));
 		if (result >= modulus)
 		{
 			result -= modulus;