using System;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks;
///
/// 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)
///
[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;
}
}
///
/// Baseline Montgomery reduction without precomputed cycle data.
///
[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;
}
///
/// Montgomery reduction with a known cycle length, excluding the cost of calculating the cycle.
///
[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;
}
///
/// Montgomery reduction with a cycle length computed on the fly using the GPU-friendly calculator.
///
[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);
}
}