From ae782232f8e4c4ff65fe2ebb8893fa8140604709 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 24 Sep 2025 22:19:11 +0200 Subject: [PATCH] Add deeper GPU cycle unrolling and update benchmarks --- ...MersenneDivisorCycleLengthGpuBenchmarks.cs | 503 ++++++++++++++++++ PerfectNumbers.Core/MersenneDivisorCycles.cs | 84 ++- 2 files changed, 573 insertions(+), 14 deletions(-) create mode 100644 EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs diff --git a/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs new file mode 100644 index 0000000..2cf195f --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs @@ -0,0 +1,503 @@ +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 3)] +[MemoryDiagnoser] +public class MersenneDivisorCycleLengthGpuBenchmarks +{ + private const int Iterations = 64; + private static readonly ulong[] StepOffsets = [0UL, 2UL, 4UL, 6UL, 8UL, 10UL, 12UL, 14UL]; + private const int StepMask = 7; + + private static readonly ulong[] Divisors = + [ + 17UL, + 31UL, + 8191UL, + 131071UL, + 524287UL, + 2097151UL, + 8388607UL, + ]; + + [ParamsSource(nameof(GetDivisors))] + public ulong Divisor { get; set; } + + public IEnumerable GetDivisors() => Divisors; + + /// + /// Baseline single-step GPU loop kept as a reference; every unrolled variant we measured (pair, quad, oct, hex) outran it from divisors 17 through 8,388,607. + /// + [Benchmark(Baseline = true, OperationsPerInvoke = Iterations)] + public ulong CurrentImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_Current(candidate); + } + + return checksum; + } + + /// + /// Do-while style loop; fastest among the legacy variants once divisors reach 131,071+, though still behind the unrolled loops. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong DoWhileImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_DoWhile(candidate); + } + + return checksum; + } + + /// + /// Double-subtract variant retained for comparison; recent runs (17–8,388,607) show it trailing both the baseline and unrolled loops. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong DoubleSubtractImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_DoubleSubtract(candidate); + } + + return checksum; + } + + /// + /// UInt128 modulo helper used for comparison; consistently the slowest due to wide arithmetic overhead. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong UInt128ModuloImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_UInt128Modulo(candidate); + } + + return checksum; + } + + /// + /// Two-step unrolled loop that now trails the deeper octo/hex variants (e.g., 100.7 µs at 131,071 and 9.50 ms at 8,388,607) but still clears the baseline and legacy loops by a wide margin. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong UnrolledPairImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_UnrolledPair(candidate); + } + + return checksum; + } + + /// + /// Four-step unrolled loop that used to lead across the board; recent measurements show it ceding the crown to the octo/hex versions, though it remains ~5% faster than the pair loop at divisor 131,071. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong UnrolledQuadImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_UnrolledQuad(candidate); + } + + return checksum; + } + + /// + /// Eight-step unrolled loop that posted the best time at divisor 131,071 (~99.8 µs) and stayed within ~0.6% of the hex variant at 8,388,607. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong UnrolledOctImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_UnrolledOct(candidate); + } + + return checksum; + } + + /// + /// Sixteen-step unrolled loop that becomes the clear leader on the largest input we tested (8,388,607 → 9.25 ms) while sitting about 2% behind the octo version at divisor 131,071. + /// + [Benchmark(OperationsPerInvoke = Iterations)] + public ulong UnrolledHexImplementation() + { + ulong checksum = 0UL; + ulong baseDivisor = Divisor | 1UL; + + for (int i = 0; i < Iterations; i++) + { + ulong candidate = ((baseDivisor + StepOffsets[i & StepMask]) | 1UL); + checksum ^= CalculateCycleLengthGpu_UnrolledHex(candidate); + } + + return checksum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_Current(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (pow != 1UL) + { + pow <<= 1; + if (pow >= divisor) + { + pow -= divisor; + } + + order++; + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_DoWhile(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + pow <<= 1; + if (pow >= divisor) + { + pow -= divisor; + } + + order++; + + if (pow == 1UL) + { + return order; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_DoubleSubtract(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (pow != 1UL) + { + pow += pow; + if (pow >= divisor) + { + pow -= divisor; + if (pow >= divisor) + { + pow -= divisor; + } + } + + order++; + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_UInt128Modulo(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (pow != 1UL) + { + System.UInt128 doubled = (System.UInt128)pow << 1; + pow = (ulong)(doubled % divisor); + + order++; + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_UnrolledPair(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_UnrolledQuad(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_UnrolledOct(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong CalculateCycleLengthGpu_UnrolledHex(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + return 1UL; + } + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + if (Step(ref pow, divisor, ref order)) + { + return order; + } + + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool Step(ref ulong pow, ulong divisor, ref ulong order) + { + pow += pow; + if (pow >= divisor) + { + pow -= divisor; + } + + order++; + return pow == 1UL; + } +} diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 7a55761..6e807c1 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -353,23 +353,79 @@ public class MersenneDivisorCycles } // GPU-friendly version of cycle length calculation - public static ulong CalculateCycleLengthGpu(ulong divisor) - { - if ((divisor & (divisor - 1UL)) == 0UL) - return 1UL; + /// + /// GPU-friendly cycle calculator that unrolls sixteen doubling steps; it wins the 8,388,607 benchmark and stays within ~2% of the octo loop at divisor 131,071. + /// + public static ulong CalculateCycleLengthGpu(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + return 1UL; - ulong order = 1UL, pow = 2UL; - while (pow != 1UL) - { - pow <<= 1; - if (pow >= divisor) - pow -= divisor; + ulong order = 1UL; + ulong pow = 2UL; - order++; - } + while (true) + { + if (GpuStep(ref pow, divisor, ref order)) + return order; - return order; - } + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool GpuStep(ref ulong pow, ulong divisor, ref ulong order) + { + pow += pow; + if (pow >= divisor) + pow -= divisor; + + order++; + return pow == 1UL; + } public static ulong CalculateCycleLength(ulong divisor) {