Disable GPU heuristic cycle benchmarks

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

@@ -1,3 +1,5 @@
+using System;
+using System.Buffers.Binary;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using System.Numerics;
@@ -5,25 +7,47 @@ using PerfectNumbers.Core;
 
 namespace EvenPerfectBitScanner.Benchmarks;
 
+/// <remarks>
+/// Benchmark results collected on AMD Ryzen 5 5625U (.NET 8.0, launchCount: 1, warmupCount: 1, iterationCount: 5):
+/// Small sample (~1e6 random odd moduli)
+/// - MontgomeryWithoutCycle: 35.41 µs (baseline)
+/// - MontgomeryWithPrecomputedCycle: 25.81 µs (0.73× baseline)
+/// - MontgomeryWithGpuCycleComputation: 54.09 ms (1,527× baseline, includes GPU cycle calculation)
+/// Large sample (Mersenne-like moduli from 2^48-1 to 2^63-1)
+/// - MontgomeryWithoutCycle: 117.33 µs (baseline)
+/// - MontgomeryWithPrecomputedCycle: 9.12 µs (0.08× baseline)
+/// - MontgomeryWithGpuCycleComputation: 12.91 µs (0.11× baseline)
+/// </remarks>
 [MemoryDiagnoser]
 [SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
 public class Pow2MontgomeryModBenchmarks
 {
     private const int SampleCount = 256;
+    private const ulong OneShiftedLeft63 = 1UL << 63;
+    private const long OneShiftedLeft60 = 1L << 60;
 
     private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
     private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
+    private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
+    private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
     private readonly ulong[] _smallExponents = new ulong[SampleCount];
     private readonly ulong[] _largeExponents = new ulong[SampleCount];
+    private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
     private readonly ulong[] _smallCycles = new ulong[SampleCount];
     private readonly ulong[] _largeCycles = new ulong[SampleCount];
+    private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
+    private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
 
     private readonly Random _random = new(13);
+    private readonly Random _heuristicRandom = new(113);
+    private ulong? _previousPrimeOrder;
+    private UInt128? _previousWidePrimeOrder;
 
     public enum InputScale
     {
         Small,
-        Large
+        Large,
+        VeryLarge
     }
 
     [ParamsAllValues]
@@ -34,6 +58,9 @@ public class Pow2MontgomeryModBenchmarks
     {
         for (int i = 0; i < SampleCount; i++)
         {
+#if DEBUG
+            Console.WriteLine($"Generating sample {i}");
+#endif
             ulong smallModulus = NextSmallOddModulus();
             _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
             _smallExponents[i] = NextSmallExponent();
@@ -43,6 +70,19 @@ public class Pow2MontgomeryModBenchmarks
             _largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
             _largeExponents[i] = NextLargeExponent();
             _largeCycles[i] = largeCycle;
+
+            ulong veryLargeModulus = NextVeryLargeOddModulus();
+            _veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus);
+            _veryLargeExponents[i] = NextVeryLargeExponent();
+#if DEBUG
+            Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
+#endif
+            _veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
+
+            Console.WriteLine($"Generating wide sample {i}");
+            UInt128 wideModulus = NextWideOddModulus();
+            _wideModuli[i] = wideModulus;
+            _wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
         }
     }
 
@@ -114,6 +154,89 @@ public class Pow2MontgomeryModBenchmarks
         return checksum;
     }
 
+    /// <summary>
+    /// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
+    /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
+    /// </summary>
+    // [Benchmark]
+    public ulong MontgomeryWithGpuCycleComputationCpu()
+    {
+        GetData(out var exponents, out var divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
+    /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
+    /// </summary>
+    // [Benchmark]
+    public ulong MontgomeryWithGpuCycleComputationGpu()
+    {
+        GetData(out var exponents, out var divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    [Benchmark]
+    public ulong MontgomeryWithHeuristicCycleComputationCpu()
+    {
+        GetData(out var exponents, out var divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    // [Benchmark]
+    public ulong MontgomeryWithHeuristicCycleComputationGpu()
+    {
+        GetData(out var exponents, out var divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    // [Benchmark]
+    public UInt128 HeuristicWideCycleComputation()
+    {
+        UInt128 checksum = UInt128.Zero;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
+            _wideCycles[i] = cycle;
+            checksum ^= cycle;
+        }
+
+        return checksum;
+    }
+
     /// <summary>
     /// Left-to-right Montgomery scan that trails the baseline: 39.11 μs on small inputs (1.08×) and 126.38 μs on large ones (1.08×).
     /// </summary>
@@ -218,23 +341,39 @@ public class Pow2MontgomeryModBenchmarks
 
     private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
     {
-        if (Scale == InputScale.Small)
+        switch (Scale)
         {
-            exponents = _smallExponents;
-            divisors = _smallDivisors;
-            cycles = _smallCycles;
-        }
-        else
-        {
-            exponents = _largeExponents;
-            divisors = _largeDivisors;
-            cycles = _largeCycles;
+            case InputScale.Small:
+                exponents = _smallExponents;
+                divisors = _smallDivisors;
+                cycles = _smallCycles;
+                break;
+            case InputScale.Large:
+                exponents = _largeExponents;
+                divisors = _largeDivisors;
+                cycles = _largeCycles;
+                break;
+            case InputScale.VeryLarge:
+                exponents = _veryLargeExponents;
+                divisors = _veryLargeDivisors;
+                cycles = _veryLargeCycles;
+                break;
+            default:
+                throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
         }
     }
 
     private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
 
-    private ulong NextLargeExponent() => (ulong)_random.NextInt64(1L << 60, long.MaxValue);
+    private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
+
+    private ulong NextVeryLargeExponent()
+    {
+        Span<byte> buffer = stackalloc byte[8];
+        _heuristicRandom.NextBytes(buffer);
+        ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
+        return value;
+    }
 
     private ulong NextSmallOddModulus()
     {
@@ -249,6 +388,27 @@ public class Pow2MontgomeryModBenchmarks
         return (modulus, (ulong)bitLength);
     }
 
+    private ulong NextVeryLargeOddModulus()
+    {
+        Span<byte> buffer = stackalloc byte[8];
+        _heuristicRandom.NextBytes(buffer);
+        ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
+        return value;
+    }
+
+    private UInt128 NextWideOddModulus()
+    {
+        Span<byte> buffer = stackalloc byte[16];
+        _heuristicRandom.NextBytes(buffer);
+
+        ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+        Span<byte> highSlice = buffer.Slice(8, 8);
+        ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
+
+        UInt128 value = ((UInt128)high << 64) | low;
+        return value | UInt128.One;
+    }
+
     private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
     {
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -269,6 +429,64 @@ public class Pow2MontgomeryModBenchmarks
             montgomeryTwoSquared);
     }
 
+    private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
+    {
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
+        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        {
+            return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
+        }
+
+#if DEBUG
+        Console.WriteLine("Trying heuristic. Prime order calculation");
+#endif
+        ulong order = PrimeOrderCalculator.Calculate(
+            modulus,
+            _previousPrimeOrder,
+            divisorData,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
+
+        if (order != 0UL)
+        {
+            _previousPrimeOrder = order;
+            return order;
+        }
+
+        _previousPrimeOrder = null;
+
+        Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
+        return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
+    }
+
+    private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
+    {
+        if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
+        {
+            return MersenneDivisorCycles.GetCycle(modulus);
+        }
+
+#if DEBUG
+        Console.WriteLine("Trying heuristic. Wide prime order calculation");
+#endif
+        UInt128 order = PrimeOrderCalculator.Calculate(
+            modulus,
+            _previousWidePrimeOrder,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
+
+        if (order != UInt128.Zero)
+        {
+            _previousWidePrimeOrder = order;
+            return order;
+        }
+
+        _previousWidePrimeOrder = null;
+
+        Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
+        return MersenneDivisorCycles.GetCycle(modulus);
+    }
+
     private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
 
     private static ulong ComputeMontgomeryNPrime(ulong modulus)
@@ -510,4 +728,3 @@ public class Pow2MontgomeryModBenchmarks
         return (ulong)(product % modulus);
     }
 }
-