From 42447c4e43050df1ddc4433f1b3676e61691ef16 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Fri, 14 Nov 2025 11:43:16 +0100 Subject: [PATCH] BinaryGcd switched to the best performing variant according to benchmarks. --- .../Gpu/GpuKernelCompilationTests.cs | 7 +- .../Gpu/Accelerators/AcceleratorStreamPool.cs | 1 - .../PrimeOrderCalculatorAccelerator.cs | 10 +- .../Gpu/Kernels/PrimeOrderKernels.cs | 18 -- .../Gpu/Kernels/PrimeTesterKernels.cs | 40 ++-- .../Gpu/PrimeOrderGpuHeuristics.cs | 175 +++++++----------- .../PrimeOrderCalculator.Cpu.cs | 2 +- PerfectNumbers.Core/PrimeTester.cs | 3 +- PerfectNumbers.Core/ULongExtensions.cs | 22 ++- 9 files changed, 112 insertions(+), 166 deletions(-) diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index d6c38f2..d3bdf9c 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -114,10 +114,6 @@ public class GpuKernelCompilationTests typeof(PrimeOrderGpuHeuristics).GetMethod("GetPow2ModKernel", BindingFlags.NonPublic | BindingFlags.Static) ?? throw new InvalidOperationException("PrimeOrderGpuHeuristics.GetPow2ModKernel not found."); - private static readonly MethodInfo PrimeOrderGetPow2ModWideKernelMethod = - typeof(PrimeOrderGpuHeuristics).GetMethod("GetPow2ModWideKernel", BindingFlags.NonPublic | BindingFlags.Static) - ?? throw new InvalidOperationException("PrimeOrderGpuHeuristics.GetPow2ModWideKernel not found."); - private static readonly Action, ArrayView1D> MersenneDivisorCyclesKernel = DivisorCycleKernels.GpuDivisorCycleKernel; @@ -325,7 +321,8 @@ public class GpuKernelCompilationTests private static void CompilePrimeOrderPow2ModWideKernel(Accelerator accelerator) { - _ = PrimeOrderGetPow2ModWideKernelMethod.Invoke(null, new object[] { accelerator }); + var kernel = accelerator.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(PrimeOrderGpuHeuristics.Pow2ModKernelWide); + _ = kernel; } private static void CompileMersenneDivisorCyclesKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs index 1a334b6..e06b040 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs +++ b/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; diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index 797faf1..2cb19b6 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/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(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, ArrayView1D, ArrayView1D>(Pow2MontgomeryKernels.TryPollardRhoKernel)); - SharesFactorKernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)); ; + Pow2ModWideKernel = KernelUtil.GetKernel( + accelerator.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(PrimeOrderGpuHeuristics.Pow2ModKernelWide)); + + SharesFactorKernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)); } public void Dispose() @@ -226,4 +232,4 @@ public sealed class PrimeOrderCalculatorAccelerator // Accelerator.Dispose(); // Context.Dispose(); } -} \ No newline at end of file +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs index 17ab513..9ac135f 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs +++ b/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, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) - { - var pool = _pow2ModWideKernel ??= []; - if (pool.TryGetValue(accelerator, out var cached)) - { - return cached; - } - - var loaded = accelerator.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(Pow2ModKernelWide); - - var kernel = KernelUtil.GetKernel(loaded); - - cached = kernel.CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); - 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 exponents, GpuUInt128 modulus, ArrayView1D remainders) { diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index a3eabff..efe1b82 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/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); + } + } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 101396e..e113f8f 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/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, MontgomeryDivisorData, ArrayView1D>>? _pow2ModKernel; - [ThreadStatic] - private static Dictionary, GpuUInt128, ArrayView1D>>? _pow2ModWideKernel; - [ThreadStatic] private static Dictionary, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>? _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(primes.Length); squareBuffer = accelerator.Allocate1D(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 primeTargets, Span 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 primeTargets, Span 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? factorBuffer; - MemoryBuffer1D? exponentBuffer; - MemoryBuffer1D? countBuffer; - MemoryBuffer1D? remainingBuffer; - MemoryBuffer1D? fullyFactoredBuffer; - // lock(accelerator) - { - factorBuffer = accelerator.Allocate1D(primeTargets.Length); - exponentBuffer = accelerator.Allocate1D(exponentTargets.Length); - countBuffer = accelerator.Allocate1D(1); - remainingBuffer = accelerator.Allocate1D(1); - fullyFactoredBuffer = accelerator.Allocate1D(1); - } + MemoryBuffer1D? factorBuffer = accelerator.Allocate1D(primeTargets.Length); + MemoryBuffer1D? exponentBuffer = accelerator.Allocate1D(exponentTargets.Length); + MemoryBuffer1D? countBuffer = accelerator.Allocate1D(1); + MemoryBuffer1D? remainingBuffer = accelerator.Allocate1D(1); + MemoryBuffer1D? fullyFactoredBuffer = accelerator.Allocate1D(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(HeuristicStackCapacity); stackProductBuffer = accelerator.Allocate1D(HeuristicStackCapacity); resultBuffer = accelerator.Allocate1D(1); - statusBuffer = accelerator.Allocate1D(1); + statusBuffer = accelerator.Allocate1D(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 exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) @@ -585,7 +563,7 @@ internal static partial class PrimeOrderGpuHeuristics // lock(accelerator) { exponentBuffer = accelerator.Allocate1D(exponents.Length); - remainderBuffer = accelerator.Allocate1D(exponents.Length); + remainderBuffer = accelerator.Allocate1D(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 exponents, UInt128 prime, Span results) + private static void ComputeOnGpuWide(ReadOnlySpan exponents, UInt128 prime, Span results) { int length = exponents.Length; - if (length == 0) - { - return true; - } GpuUInt128[]? rentedExponents = null; GpuUInt128[]? rentedResults = null; - // GpuPrimeWorkLimiter.Acquire(); ArrayPool gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool; + PrimeOrderCalculatorAccelerator? gpu = null; + bool poolingRequired = length > WideStackThreshold; + Span exponentSpan = poolingRequired + ? stackalloc GpuUInt128[length] + : new Span(rentedExponents = gpuUInt128Pool.Rent(length), 0, length); + Span resultSpan = length <= WideStackThreshold + ? stackalloc GpuUInt128[length] + : new Span(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? exponentBuffer; - MemoryBuffer1D? remainderBuffer; - // lock(accelerator) - { - exponentBuffer = accelerator.Allocate1D(length); - remainderBuffer = accelerator.Allocate1D(length); - } - - Span exponentSpan = length <= WideStackThreshold - ? stackalloc GpuUInt128[length] - : new Span(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 resultSpan = length <= WideStackThreshold - ? stackalloc GpuUInt128[length] - : new Span(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 exponentBuffer = accelerator.Allocate1D(length); + MemoryBuffer1D remainderBuffer = accelerator.Allocate1D(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); } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 2b4aea0..597e69e 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/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 pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index c357f65..bc5ccd5 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/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(); diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index a0726e8..bfdcb3e 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/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.