Merge pull request #94 from mlemanczyk/codex/check-current-benchmarks-for-montgomery-method

Document cycle computation benchmark results

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

@@ -0,0 +1,174 @@
+using System;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Jobs;
+using PerfectNumbers.Core;
+
+namespace EvenPerfectBitScanner.Benchmarks;
+
+/// <remarks>
+/// Benchmark results collected on Intel Xeon Platinum 8370C (Ubuntu 24.04, .NET 8.0, launchCount: 1, warmupCount: 1,
+/// iterationCount: 5):
+/// Small sample (~1e6 exponents and random odd moduli)
+/// - MontgomeryWithoutCycle: 43.17 µs (baseline)
+/// - MontgomeryWithPrecomputedCycle: 29.88 µs (0.69× baseline)
+/// - MontgomeryWithGpuCycleComputation: 75.09 ms (1740× baseline, includes GPU cycle calculation)
+/// Large sample (Mersenne-like moduli from 2^48-1 to 2^63-1)
+/// - MontgomeryWithoutCycle: 145.32 µs (baseline)
+/// - MontgomeryWithPrecomputedCycle: 10.93 µs (0.08× baseline)
+/// - MontgomeryWithGpuCycleComputation: 16.55 µs (0.11× baseline, includes GPU cycle calculation)
+/// </remarks>
+[MemoryDiagnoser]
+[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
+public class Pow2MontgomeryModCycleComputationBenchmarks
+{
+    private const int SampleCount = 256;
+
+    private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
+    private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
+    private readonly ulong[] _smallExponents = new ulong[SampleCount];
+    private readonly ulong[] _largeExponents = new ulong[SampleCount];
+    private readonly ulong[] _smallCycles = new ulong[SampleCount];
+    private readonly ulong[] _largeCycles = new ulong[SampleCount];
+
+    private readonly Random _random = new(113);
+
+    public enum InputScale
+    {
+        Small,
+        Large
+    }
+
+    [ParamsAllValues]
+    public InputScale Scale { get; set; }
+
+    [GlobalSetup]
+    public void Setup()
+    {
+        for (int i = 0; i < SampleCount; i++)
+        {
+            ulong smallModulus = NextSmallOddModulus();
+            _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
+            _smallExponents[i] = NextSmallExponent();
+            _smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus);
+
+            (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
+            _largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
+            _largeExponents[i] = NextLargeExponent();
+            _largeCycles[i] = largeCycle;
+        }
+    }
+
+    /// <summary>
+    /// Baseline Montgomery reduction without precomputed cycle data.
+    /// </summary>
+    [Benchmark(Baseline = true)]
+    public ulong MontgomeryWithoutCycle()
+    {
+        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            checksum ^= exponents[i].Pow2MontgomeryMod(divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// Montgomery reduction with a known cycle length, excluding the cost of calculating the cycle.
+    /// </summary>
+    [Benchmark]
+    public ulong MontgomeryWithPrecomputedCycle()
+    {
+        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// Montgomery reduction with a cycle length computed on the fly using the GPU-friendly calculator.
+    /// </summary>
+    [Benchmark]
+    public ulong MontgomeryWithGpuCycleComputation()
+    {
+        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
+        }
+
+        return checksum;
+    }
+
+    private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
+    {
+        if (Scale == InputScale.Small)
+        {
+            exponents = _smallExponents;
+            divisors = _smallDivisors;
+            cycles = _smallCycles;
+        }
+        else
+        {
+            exponents = _largeExponents;
+            divisors = _largeDivisors;
+            cycles = _largeCycles;
+        }
+    }
+
+    private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
+
+    private ulong NextLargeExponent() => (ulong)_random.NextInt64(1L << 60, long.MaxValue);
+
+    private ulong NextSmallOddModulus()
+    {
+        ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
+        return value | 1UL;
+    }
+
+    private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
+    {
+        int bitLength = _random.Next(48, 64);
+        ulong modulus = (1UL << bitLength) - 1UL;
+        return (modulus, (ulong)bitLength);
+    }
+
+    private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
+    {
+        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        {
+            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
+        }
+
+        return new MontgomeryDivisorData(
+            modulus,
+            ComputeMontgomeryNPrime(modulus),
+            ComputeMontgomeryResidue(1UL, modulus),
+            ComputeMontgomeryResidue(2UL, modulus));
+    }
+
+    private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
+
+    private static ulong ComputeMontgomeryNPrime(ulong modulus)
+    {
+        ulong inv = modulus;
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        return unchecked(0UL - inv);
+    }
+}
+