BinaryGcd switched to the best performing variant according to benchmarks.

PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs

@@ -1,6 +1,5 @@
 using System.Collections.Concurrent;
 using ILGPU.Runtime;
-using Open.Collections;
 
 namespace PerfectNumbers.Core.Gpu.Accelerators;
 

PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs

@@ -1,6 +1,7 @@
 using System.Runtime.CompilerServices;
 using ILGPU;
 using ILGPU.Runtime;
+using PerfectNumbers.Core.Gpu;
 using PerfectNumbers.Core.Gpu.Kernels;
 
 namespace PerfectNumbers.Core.Gpu.Accelerators;
@@ -110,6 +111,7 @@ public sealed class PrimeOrderCalculatorAccelerator
 	public readonly Kernel ConvertToStandardKernel;
 	public readonly Kernel KeepMontgomeryKernel;
 	public readonly Kernel PollardRhoKernel;
+	public readonly Kernel Pow2ModWideKernel;
 	public Kernel SharesFactorKernel;
 
 
@@ -180,6 +182,7 @@ public sealed class PrimeOrderCalculatorAccelerator
 		SpecialMaxResult = accelerator.Allocate1D<ulong>(1);
 
 		var sharedTables = LastDigitGpuTables.EnsureStaticTables(accelerator);
+
 		DevicePrimesLastOne = sharedTables.DevicePrimesLastOne;
 		DevicePrimesLastSeven = sharedTables.DevicePrimesLastSeven;
 		DevicePrimesLastThree = sharedTables.DevicePrimesLastThree;
@@ -202,7 +205,10 @@ public sealed class PrimeOrderCalculatorAccelerator
 
 		PollardRhoKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel<ulong, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<byte, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2MontgomeryKernels.TryPollardRhoKernel));
 
-		SharesFactorKernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SharesFactorKernel)); ;
+		Pow2ModWideKernel = KernelUtil.GetKernel(
+			accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>>(PrimeOrderGpuHeuristics.Pow2ModKernelWide));
+
+		SharesFactorKernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SharesFactorKernel));
 	}
 
 	public void Dispose()
@@ -226,4 +232,4 @@ public sealed class PrimeOrderCalculatorAccelerator
 		// Accelerator.Dispose();
 		// Context.Dispose();
 	}
-}
+}

PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs

@@ -2,7 +2,6 @@ using System.Numerics;
 using System.Runtime.CompilerServices;
 using ILGPU;
 using ILGPU.Runtime;
-using PerfectNumbers.Core.Cpu;
 
 namespace PerfectNumbers.Core.Gpu;
 
@@ -1207,23 +1206,6 @@ internal static partial class PrimeOrderGpuHeuristics
 		remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent);
 	}
 
-	private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)
-	{
-		var pool = _pow2ModWideKernel ??= [];
-		if (pool.TryGetValue(accelerator, out var cached))
-		{
-			return cached;
-		}
-
-		var loaded = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>>(Pow2ModKernelWide);
-
-		var kernel = KernelUtil.GetKernel(loaded);
-
-		cached = kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>>>();
-		pool[accelerator] = cached;
-		return cached;
-	}
-
 	/// This kernel always sets the result of the corresponding element. Callers don't need to clear the output buffers.
 	internal static void Pow2ModKernelWide(Index1D index, ArrayView1D<GpuUInt128, Stride1D.Dense> exponents, GpuUInt128 modulus, ArrayView1D<GpuUInt128, Stride1D.Dense> remainders)
 	{

PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs

@@ -121,35 +121,37 @@ internal static class PrimeTesterKernels
 
     private static ulong BinaryGcdGpu(ulong u, ulong v)
     {
-        // TODO: Replace this inline GPU binary GCD with the kernel extracted from
-        // GpuUInt128BinaryGcdBenchmarks via GpuKernelPool so device callers reuse the
-        // fully unrolled ladder instead of this branchy fallback.
         bool uIsZero = u == 0UL;
         bool vIsZero = v == 0UL;
         bool eitherZero = uIsZero | vIsZero;
         ulong zeroResult = uIsZero ? v : u;
 
-        ulong combined = u | v | 1UL;
+        ulong combined = (u | v) | 1UL;
         int shift = XMath.TrailingZeroCount(combined);
+        ulong currentU = u >> XMath.TrailingZeroCount(u | 1UL);
+        ulong currentV = v >> XMath.TrailingZeroCount(v | 1UL);
 
-        ulong normalizedU = u >> XMath.TrailingZeroCount(u | 1UL);
-        ulong normalizedV = v >> XMath.TrailingZeroCount(v | 1UL);
-
-        ulong currentU = normalizedU;
-        ulong currentV = normalizedV;
-        bool loopCondition = currentV != 0UL;
-
-        while (loopCondition)
+        while (true)
         {
             currentV >>= XMath.TrailingZeroCount(currentV);
-            ulong minValue = XMath.Min(currentU, currentV);
-            ulong maxValue = XMath.Max(currentU, currentV);
-            currentU = minValue;
+            bool swap = currentU > currentV;
+            ulong minValue = ConditionalSelect(currentU, currentV, swap);
+            ulong maxValue = ConditionalSelect(currentV, currentU, swap);
             currentV = maxValue - minValue;
-            loopCondition = currentV != 0UL;
+            currentU = minValue;
+            if (currentV == 0UL)
+            {
+                ulong gcd = currentU << shift;
+                return eitherZero ? zeroResult : gcd;
+            }
         }
-
-        ulong gcd = currentU << shift;
-        return eitherZero ? zeroResult : gcd;
     }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong ConditionalSelect(ulong left, ulong right, bool useRight)
+    {
+        ulong mask = useRight ? ulong.MaxValue : 0UL;
+        return left ^ ((left ^ right) & mask);
+    }
+
 }

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -1,9 +1,10 @@
 using System.Buffers;
 using System.Collections.Concurrent;
+using System.Runtime.InteropServices;
 using ILGPU;
 using ILGPU.Runtime;
-using System.Runtime.InteropServices;
 using ILGPU.Runtime.OpenCL;
+using PerfectNumbers.Core.Gpu.Accelerators;
 namespace PerfectNumbers.Core.Gpu;
 
 internal enum GpuPow2ModStatus
@@ -26,9 +27,6 @@ internal static partial class PrimeOrderGpuHeuristics
 	[ThreadStatic]
 	private static Dictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>>>? _pow2ModKernel;
 
-	[ThreadStatic]
-	private static Dictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>>>? _pow2ModWideKernel;
-
 	[ThreadStatic]
 	private static Dictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, int, ulong, uint, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<byte, Stride1D.Dense>>>? _partialFactorKernel;
 
@@ -89,7 +87,7 @@ internal static partial class PrimeOrderGpuHeuristics
 		public int Count;
 	}
 
-	private const int WideStackThreshold = 8;
+	private const int WideStackThreshold = 12;
 	private static PrimeOrderGpuCapability s_capability = PrimeOrderGpuCapability.Default;
 
 	private const int Pow2WindowSizeBits = 8;
@@ -130,7 +128,7 @@ internal static partial class PrimeOrderGpuHeuristics
 			primeBuffer = accelerator.Allocate1D<uint>(primes.Length);
 			squareBuffer = accelerator.Allocate1D<ulong>(squares.Length);
 		}
-		
+
 		primeBuffer.View.CopyFromCPU(stream, primes);
 		squareBuffer.View.CopyFromCPU(stream, squares);
 		cached = new SmallPrimeDeviceCache
@@ -163,17 +161,16 @@ internal static partial class PrimeOrderGpuHeuristics
 	}
 
 	public static bool TryPartialFactor(
+		PrimeOrderCalculatorAccelerator gpu,
 		ulong value,
 		uint limit,
 		Span<ulong> primeTargets,
 		Span<int> exponentTargets,
 		out int factorCount,
-		out ulong remaining,
-		out bool fullyFactored)
+		out ulong remaining)
 	{
 		factorCount = 0;
 		remaining = value;
-		fullyFactored = false;
 
 		if (primeTargets.Length == 0 || exponentTargets.Length == 0)
 		{
@@ -184,13 +181,13 @@ internal static partial class PrimeOrderGpuHeuristics
 		exponentTargets.Clear();
 
 		if (!TryLaunchPartialFactorKernel(
+				gpu,
 				value,
 				limit,
 				primeTargets,
 				exponentTargets,
 				out int extracted,
-				out ulong leftover,
-				out bool kernelFullyFactored))
+				out ulong leftover))
 		{
 			return false;
 		}
@@ -200,28 +197,25 @@ internal static partial class PrimeOrderGpuHeuristics
 		{
 			factorCount = 0;
 			remaining = value;
-			fullyFactored = false;
 			return false;
 		}
 
 		factorCount = extracted;
 		remaining = leftover;
-		fullyFactored = kernelFullyFactored && leftover == 1UL;
 		return true;
 	}
 
 	private static bool TryLaunchPartialFactorKernel(
+		PrimeOrderCalculatorAccelerator gpu,
 		ulong value,
 		uint limit,
 		Span<ulong> primeTargets,
 		Span<int> exponentTargets,
 		out int factorCount,
-		out ulong remaining,
-		out bool fullyFactored)
+		out ulong remaining)
 	{
 		factorCount = 0;
 		remaining = value;
-		fullyFactored = false;
 
 		try
 		{
@@ -232,21 +226,13 @@ internal static partial class PrimeOrderGpuHeuristics
 			SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator, stream);
 
 			// TODO: We should create / reallocate these buffer only once or if the new required length is bigger than capacity.
-			MemoryBuffer1D<ulong, Stride1D.Dense>? factorBuffer;
-			MemoryBuffer1D<int, Stride1D.Dense>? exponentBuffer;
-			MemoryBuffer1D<int, Stride1D.Dense>? countBuffer;
-			MemoryBuffer1D<ulong, Stride1D.Dense>? remainingBuffer;
-			MemoryBuffer1D<byte, Stride1D.Dense>? fullyFactoredBuffer;
-			// lock(accelerator)
-			{
-				factorBuffer = accelerator.Allocate1D<ulong>(primeTargets.Length);
-				exponentBuffer = accelerator.Allocate1D<int>(exponentTargets.Length);
-				countBuffer = accelerator.Allocate1D<int>(1);
-				remainingBuffer = accelerator.Allocate1D<ulong>(1);
-				fullyFactoredBuffer = accelerator.Allocate1D<byte>(1);				
-			}
+			MemoryBuffer1D<ulong, Stride1D.Dense>? factorBuffer = accelerator.Allocate1D<ulong>(primeTargets.Length);
+			MemoryBuffer1D<int, Stride1D.Dense>? exponentBuffer = accelerator.Allocate1D<int>(exponentTargets.Length);
+			MemoryBuffer1D<int, Stride1D.Dense>? countBuffer = accelerator.Allocate1D<int>(1);
+			MemoryBuffer1D<ulong, Stride1D.Dense>? remainingBuffer = accelerator.Allocate1D<ulong>(1);
+			MemoryBuffer1D<byte, Stride1D.Dense>? fullyFactoredBuffer = accelerator.Allocate1D<byte>(1);
 
-			// TODO: There is no need to clear these buffers because the kernel will always assign values within the required bounds.
+			// There is no need to clear these buffers because the kernel will always assign values within the required bounds. Keep it commented out.
 			// factorBuffer.MemSetToZero(stream);
 			// exponentBuffer.MemSetToZero(stream);
 			// countBuffer.MemSetToZero(stream);
@@ -281,7 +267,6 @@ internal static partial class PrimeOrderGpuHeuristics
 			stream.Synchronize();
 			stream.Dispose();
 
-			fullyFactored = fullyFactoredFlag != 0;
 			factorCount = Math.Min(factorCount, primeTargets.Length);
 
 			// TODO: These buffers shouldn't be disposed but rather reused accross call with the assigned accelerator.
@@ -296,18 +281,12 @@ internal static partial class PrimeOrderGpuHeuristics
 		catch (CLException ex)
 		{
 			Console.WriteLine($"GPU ERROR ({ex.Error}): {ex.Message}");
-			factorCount = 0;
-			remaining = value;
-			fullyFactored = false;
-			return false;
+			throw;
 		}
 		catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException)
 		{
 			Console.WriteLine($"GPU ERROR: {ex.Message}");
-			factorCount = 0;
-			remaining = value;
-			fullyFactored = false;
-			return false;
+			throw;
 		}
 	}
 
@@ -404,8 +383,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
 		// GpuPrimeWorkLimiter.Acquire();
 		Accelerator accelerator = AcceleratorPool.Shared.Rent();
-		AcceleratorStream stream = accelerator.CreateStream();
-
+		AcceleratorStream stream = AcceleratorStreamPool.Rent(accelerator);
 		SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator, stream);
 
 		// TODO: These buffers should be created reallocated once and assigned to an accelerator. Callers should use their own
@@ -433,7 +411,7 @@ internal static partial class PrimeOrderGpuHeuristics
 			stackExponentBuffer = accelerator.Allocate1D<int>(HeuristicStackCapacity);
 			stackProductBuffer = accelerator.Allocate1D<ulong>(HeuristicStackCapacity);
 			resultBuffer = accelerator.Allocate1D<ulong>(1);
-			statusBuffer = accelerator.Allocate1D<byte>(1);			
+			statusBuffer = accelerator.Allocate1D<byte>(1);
 		}
 
 		// TODO: Remove the cleaning after the order kernel is modified to always set the result.
@@ -568,8 +546,8 @@ internal static partial class PrimeOrderGpuHeuristics
 			}
 		}
 
-		bool computed = TryComputeOnGpuWide(exponents, prime, target);
-		return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable;
+		ComputeOnGpuWide(exponents, prime, target);
+		return GpuPow2ModStatus.Success;
 	}
 
 	private static bool TryComputeOnGpu(ReadOnlySpan<ulong> exponents, ulong prime, in MontgomeryDivisorData divisorData, Span<ulong> results)
@@ -585,7 +563,7 @@ internal static partial class PrimeOrderGpuHeuristics
 		// lock(accelerator)
 		{
 			exponentBuffer = accelerator.Allocate1D<ulong>(exponents.Length);
-			remainderBuffer = accelerator.Allocate1D<ulong>(exponents.Length);			
+			remainderBuffer = accelerator.Allocate1D<ulong>(exponents.Length);
 		}
 
 		exponentBuffer.View.CopyFromCPU(stream, ref MemoryMarshal.GetReference(exponents), exponents.Length);
@@ -628,83 +606,56 @@ internal static partial class PrimeOrderGpuHeuristics
 		}
 	}
 
-	private static bool TryComputeOnGpuWide(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> results)
+	private static void ComputeOnGpuWide(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> results)
 	{
 		int length = exponents.Length;
-		if (length == 0)
-		{
-			return true;
-		}
 
 		GpuUInt128[]? rentedExponents = null;
 		GpuUInt128[]? rentedResults = null;
-		// GpuPrimeWorkLimiter.Acquire();
 		ArrayPool<GpuUInt128> gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool;
+		PrimeOrderCalculatorAccelerator? gpu = null;
+		bool poolingRequired = length > WideStackThreshold;
+		Span<GpuUInt128> exponentSpan = poolingRequired
+			? stackalloc GpuUInt128[length]
+			: new Span<GpuUInt128>(rentedExponents = gpuUInt128Pool.Rent(length), 0, length);
+		Span<GpuUInt128> resultSpan = length <= WideStackThreshold
+			? stackalloc GpuUInt128[length]
+			: new Span<GpuUInt128>(rentedResults = gpuUInt128Pool.Rent(length), 0, length);
 
-		try
+
+		for (int i = 0; i < length; i++)
 		{
-			Accelerator accelerator = AcceleratorPool.Shared.Rent();
-			AcceleratorStream stream = accelerator.CreateStream();
-			// TODO: These buffers should be allocated once per accelerator and only reallocated when the new length exceeds capacity.
-			// Modify the callers to use their own pool of buffers per accelerator, so that other threads don't use the accelerators
-			// with out preallocated buffers. Share the pool with Pow2Mod kernel.
-			MemoryBuffer1D<GpuUInt128, Stride1D.Dense>? exponentBuffer;
-			MemoryBuffer1D<GpuUInt128, Stride1D.Dense>? remainderBuffer;
-			// lock(accelerator)
-			{
-				exponentBuffer = accelerator.Allocate1D<GpuUInt128>(length);
-				remainderBuffer = accelerator.Allocate1D<GpuUInt128>(length);				
-			}
-
-			Span<GpuUInt128> exponentSpan = length <= WideStackThreshold
-				? stackalloc GpuUInt128[length]
-				: new Span<GpuUInt128>(rentedExponents = gpuUInt128Pool.Rent(length), 0, length);
-
-			for (int i = 0; i < length; i++)
-			{
-				exponentSpan[i] = (GpuUInt128)exponents[i];
-			}
-
-			exponentBuffer.View.CopyFromCPU(stream, ref MemoryMarshal.GetReference(exponentSpan), length);
-			// Pow2ModWide kernel always assigns the required output elements so we don't need to worry about clearing these.
-			// remainderBuffer.MemSetToZero(stream);
-
-			GpuUInt128 modulus = (GpuUInt128)prime;
-			var pow2ModWideKernel = GetPow2ModWideKernel(accelerator);
-			pow2ModWideKernel(stream, length, exponentBuffer.View, modulus, remainderBuffer.View);
-
-
-			Span<GpuUInt128> resultSpan = length <= WideStackThreshold
-				? stackalloc GpuUInt128[length]
-				: new Span<GpuUInt128>(rentedResults = gpuUInt128Pool.Rent(length), 0, length);
-
-			remainderBuffer.View.CopyToCPU(stream, ref MemoryMarshal.GetReference(resultSpan), length);
-			stream.Synchronize();
-			stream.Dispose();
-
-			for (int i = 0; i < length; i++)
-			{
-				results[i] = (UInt128)resultSpan[i];
-			}
-
-			exponentBuffer.Dispose();
-			remainderBuffer.Dispose();
-			return true;
+			exponentSpan[i] = (GpuUInt128)exponents[i];
 		}
-		catch (Exception)
-		{
-			Console.WriteLine($"Exception when computing {prime} with exponents: {string.Join(", ", exponents.ToArray())}");
-			throw;
-		}
-		finally
-		{
-			if (rentedExponents is not null)
-			{
-				gpuUInt128Pool.Return(rentedExponents, clearArray: false);
-				gpuUInt128Pool.Return(rentedResults!, clearArray: false);
-			}
 
-			// GpuPrimeWorkLimiter.Release();
+		GpuUInt128 modulus = (GpuUInt128)prime;
+		gpu = PrimeOrderCalculatorAccelerator.Rent(1);
+		Accelerator accelerator = gpu.Accelerator;
+
+		MemoryBuffer1D<GpuUInt128, Stride1D.Dense> exponentBuffer = accelerator.Allocate1D<GpuUInt128>(length);
+		MemoryBuffer1D<GpuUInt128, Stride1D.Dense> remainderBuffer = accelerator.Allocate1D<GpuUInt128>(length);
+
+		AcceleratorStream? stream = AcceleratorStreamPool.Rent(accelerator);
+		exponentBuffer.View.CopyFromCPU(stream, exponentSpan);
+
+		gpu.Pow2ModWideKernel.Launch(stream, length, exponentBuffer.View, modulus, remainderBuffer.View);
+		remainderBuffer.View.CopyToCPU(stream, ref MemoryMarshal.GetReference(resultSpan), length);
+		stream.Synchronize();
+
+		AcceleratorStreamPool.Return(stream);
+		for (int i = 0; i < length; i++)
+		{
+			results[i] = (UInt128)resultSpan[i];
+		}
+
+		remainderBuffer.Dispose();
+		exponentBuffer.Dispose();
+		PrimeOrderCalculatorAccelerator.Return(gpu);
+
+		if (rentedExponents is not null)
+		{
+			gpuUInt128Pool.Return(rentedExponents, clearArray: false);
+			gpuUInt128Pool.Return(rentedResults!, clearArray: false);
 		}
 	}
 

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -1109,7 +1109,7 @@ internal static partial class PrimeOrderCalculator
 		bool gpuFactored = false;
 		if (IsGpuHeuristicDevice)
 		{
-			gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
+			gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out factorCount, out remaining);
 		}
 
 		List<PendingEntry> pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2);

PerfectNumbers.Core/PrimeTester.cs

@@ -289,7 +289,8 @@ public sealed class PrimeTester
 		AcceleratorStream stream = AcceleratorStreamPool.Rent(accelerator);
 		inputBuffer.View.CopyFromCPU(stream, values);
 
-		kernel.Launch(stream, 1, length, inputBuffer.View, resultBuffer.View);
+		kernel.Launch(stream, length, inputBuffer.View, resultBuffer.View);
+
 		resultBuffer.View.CopyToCPU(stream, in results);
 		stream.Synchronize();
 

PerfectNumbers.Core/ULongExtensions.cs

@@ -314,20 +314,28 @@ public static partial class ULongExtensions
 		int shift = BitOperations.TrailingZeroCount(a | b);
 		a >>= BitOperations.TrailingZeroCount(a);
 
-		do
+		while (true)
 		{
 			b >>= BitOperations.TrailingZeroCount(b);
 
-			if (a > b)
+			bool swap = a > b;
+			ulong min = Select(a, b, swap);
+			ulong max = Select(b, a, swap);
+			b = max - min;
+			a = min;
+
+			if (b == 0UL)
 			{
-				(a, b) = (b, a);
+				return a << shift;
 			}
-
-			b -= a;
 		}
-		while (b != 0UL);
+	}
 
-		return a << shift;
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static ulong Select(ulong left, ulong right, bool useRight)
+	{
+		ulong mask = useRight ? ulong.MaxValue : 0UL;
+		return left ^ ((left ^ right) & mask);
 	}
 
 	// Benchmarks (Mod5ULongBenchmarks) show the direct `% 5` is still cheaper (~0.26 ns vs 0.43 ns), so keep the modulo until a faster lookup is proven.