Formatting fix

This commit is contained in:
mlemanczyk committed 2025-10-10 03:06:01 +02:00
1 parent 0be9596bd5
commit 70eca5126b
1 file changed
+259 -259
@@ -22,321 +22,321 @@ namespace EvenPerfectBitScanner.Benchmarks;
[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] [SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
public class Pow2MontgomeryModCycleComputationBenchmarks public class Pow2MontgomeryModCycleComputationBenchmarks
{ {
private const int SampleCount = 256; private const int SampleCount = 256;
private const ulong OneShiftedLeft63 = 1UL << 63; private const ulong OneShiftedLeft63 = 1UL << 63;
private const long OneShiftedLeft60 = 1L << 60; private const long OneShiftedLeft60 = 1L << 60;
private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount]; private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount]; private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount]; private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
private readonly UInt128[] _wideModuli = new UInt128[SampleCount]; private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
private readonly ulong[] _smallExponents = new ulong[SampleCount]; private readonly ulong[] _smallExponents = new ulong[SampleCount];
private readonly ulong[] _largeExponents = new ulong[SampleCount]; private readonly ulong[] _largeExponents = new ulong[SampleCount];
private readonly ulong[] _veryLargeExponents = new ulong[SampleCount]; private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
private readonly ulong[] _smallCycles = new ulong[SampleCount]; private readonly ulong[] _smallCycles = new ulong[SampleCount];
private readonly ulong[] _largeCycles = new ulong[SampleCount]; private readonly ulong[] _largeCycles = new ulong[SampleCount];
private readonly ulong[] _veryLargeCycles = new ulong[SampleCount]; private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
private readonly UInt128[] _wideCycles = new UInt128[SampleCount]; private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
private readonly Random _random = new(113); private readonly Random _random = new(113);
private ulong? _previousPrimeOrder; private ulong? _previousPrimeOrder;
private UInt128? _previousWidePrimeOrder; private UInt128? _previousWidePrimeOrder;
public enum InputScale public enum InputScale
{ {
Small, Small,
Large, Large,
VeryLarge VeryLarge
} }
[ParamsAllValues] [ParamsAllValues]
public InputScale Scale { get; set; } public InputScale Scale { get; set; }
[GlobalSetup] [GlobalSetup]
public void Setup() public void Setup()
{ {
for (int i = 0; i < SampleCount; i++) for (int i = 0; i < SampleCount; i++)
{ {
#if DEBUG #if DEBUG
Console.WriteLine($"Generating sample {i}"); Console.WriteLine($"Generating sample {i}");
#endif #endif
ulong smallModulus = NextSmallOddModulus(); ulong smallModulus = NextSmallOddModulus();
_smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
_smallExponents[i] = NextSmallExponent(); _smallExponents[i] = NextSmallExponent();
_smallCycles[i] = CalculateCycleLengthWithHeuristics(smallModulus); _smallCycles[i] = CalculateCycleLengthWithHeuristics(smallModulus);
(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(); (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
_largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus); _largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
_largeExponents[i] = NextLargeExponent(); _largeExponents[i] = NextLargeExponent();
_largeCycles[i] = CalculateCycleLengthWithHeuristics(largeModulus, largeCycle); _largeCycles[i] = CalculateCycleLengthWithHeuristics(largeModulus, largeCycle);
ulong veryLargeModulus = NextVeryLargeOddModulus(); ulong veryLargeModulus = NextVeryLargeOddModulus();
_veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus); _veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus);
_veryLargeExponents[i] = NextVeryLargeExponent(); _veryLargeExponents[i] = NextVeryLargeExponent();
#if DEBUG #if DEBUG
Console.WriteLine($"Calculating cycle length {veryLargeModulus}"); Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
#endif #endif
_veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus); _veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
Console.WriteLine($"Generating wide sample {i}"); Console.WriteLine($"Generating wide sample {i}");
UInt128 wideModulus = NextWideOddModulus(); UInt128 wideModulus = NextWideOddModulus();
_wideModuli[i] = wideModulus; _wideModuli[i] = wideModulus;
_wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus); _wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
} }
} }
/// <summary> /// <summary>
/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
/// large set. /// large set.
/// </summary> /// </summary>
// [Benchmark(Baseline = true)] // [Benchmark(Baseline = true)]
public ulong MontgomeryWithoutCycle() public ulong MontgomeryWithoutCycle()
{ {
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL; ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++) for (int i = 0; i < SampleCount; i++)
{ {
checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false); checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
} }
return checksum; return checksum;
} }
/// <summary> /// <summary>
/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs /// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
/// on the large one (0.08×). /// on the large one (0.08×).
/// </summary> /// </summary>
[Benchmark] [Benchmark]
public ulong MontgomeryWithPrecomputedCycle() public ulong MontgomeryWithPrecomputedCycle()
{ {
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles); GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
ulong checksum = 0UL; ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++) for (int i = 0; i < SampleCount; i++)
{ {
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]); checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
} }
return checksum; return checksum;
} }
/// <summary> /// <summary>
/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to /// 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). /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
/// </summary> /// </summary>
// [Benchmark] // [Benchmark]
public ulong MontgomeryWithGpuCycleComputation() public ulong MontgomeryWithGpuCycleComputation()
{ {
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL; ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++) for (int i = 0; i < SampleCount; i++)
{ {
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus); ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]); checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
} }
return checksum; return checksum;
} }
[Benchmark] [Benchmark]
public ulong MontgomeryWithHeuristicCycleComputation() public ulong MontgomeryWithHeuristicCycleComputation()
{ {
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL; ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++) for (int i = 0; i < SampleCount; i++)
{ {
ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus); ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]); checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
} }
return checksum; return checksum;
} }
[Benchmark] [Benchmark]
public UInt128 HeuristicWideCycleComputation() public UInt128 HeuristicWideCycleComputation()
{ {
UInt128 checksum = UInt128.Zero; UInt128 checksum = UInt128.Zero;
for (int i = 0; i < SampleCount; i++) for (int i = 0; i < SampleCount; i++)
{ {
UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]); UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
_wideCycles[i] = cycle; _wideCycles[i] = cycle;
checksum ^= cycle; checksum ^= cycle;
} }
return checksum; return checksum;
} }
private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles) private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
{ {
switch (Scale) switch (Scale)
{ {
case InputScale.Small: case InputScale.Small:
exponents = _smallExponents; exponents = _smallExponents;
divisors = _smallDivisors; divisors = _smallDivisors;
cycles = _smallCycles; cycles = _smallCycles;
break; break;
case InputScale.Large: case InputScale.Large:
exponents = _largeExponents; exponents = _largeExponents;
divisors = _largeDivisors; divisors = _largeDivisors;
cycles = _largeCycles; cycles = _largeCycles;
break; break;
case InputScale.VeryLarge: case InputScale.VeryLarge:
exponents = _veryLargeExponents; exponents = _veryLargeExponents;
divisors = _veryLargeDivisors; divisors = _veryLargeDivisors;
cycles = _veryLargeCycles; cycles = _veryLargeCycles;
break; break;
default: default:
throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null); throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
} }
} }
private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L); private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue); private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
private ulong NextVeryLargeExponent() private ulong NextVeryLargeExponent()
{ {
Span<byte> buffer = stackalloc byte[8]; Span<byte> buffer = stackalloc byte[8];
_random.NextBytes(buffer); _random.NextBytes(buffer);
ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63); ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
return value; return value;
} }
private ulong NextSmallOddModulus() private ulong NextSmallOddModulus()
{ {
ulong value = (ulong)_random.NextInt64(3L, 1_000_000L); ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
return value | 1UL; return value | 1UL;
} }
private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle() private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
{ {
int bitLength = _random.Next(48, 64); int bitLength = _random.Next(48, 64);
ulong modulus = (1UL << bitLength) - 1UL; ulong modulus = (1UL << bitLength) - 1UL;
return (modulus, (ulong)bitLength); return (modulus, (ulong)bitLength);
} }
private ulong NextVeryLargeOddModulus() private ulong NextVeryLargeOddModulus()
{ {
Span<byte> buffer = stackalloc byte[8]; Span<byte> buffer = stackalloc byte[8];
_random.NextBytes(buffer); _random.NextBytes(buffer);
ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL; ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
return value; return value;
} }
private UInt128 NextWideOddModulus() private UInt128 NextWideOddModulus()
{ {
Span<byte> buffer = stackalloc byte[16]; Span<byte> buffer = stackalloc byte[16];
_random.NextBytes(buffer); _random.NextBytes(buffer);
ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
Span<byte> highSlice = buffer.Slice(8, 8); Span<byte> highSlice = buffer.Slice(8, 8);
ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63); ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
UInt128 value = ((UInt128)high << 64) | low; UInt128 value = ((UInt128)high << 64) | low;
return value | UInt128.One; return value | UInt128.One;
} }
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
{ {
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL); return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
} }
ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus); ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData( return new MontgomeryDivisorData(
modulus, modulus,
nPrime, nPrime,
montgomeryOne, montgomeryOne,
montgomeryTwo, montgomeryTwo,
montgomeryTwoSquared); montgomeryTwoSquared);
} }
private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL) private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData); return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
} }
// When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set. // When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus);
// // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None); // // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None);
// if (!isPrime) // if (!isPrime)
// { // {
// return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus); // return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
// } // }
#if DEBUG #if DEBUG
Console.WriteLine("Trying heuristic. Prime order calculation"); Console.WriteLine("Trying heuristic. Prime order calculation");
#endif #endif
ulong order = PrimeOrderCalculator.Calculate( ulong order = PrimeOrderCalculator.Calculate(
modulus, modulus,
_previousPrimeOrder, _previousPrimeOrder,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
if (order != 0UL) if (order != 0UL)
{ {
_previousPrimeOrder = order; _previousPrimeOrder = order;
return order; return order;
} }
_previousPrimeOrder = null; _previousPrimeOrder = null;
Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation"); Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData); return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
} }
private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus) private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
{ {
if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero) if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
{ {
return MersenneDivisorCycles.GetCycle(modulus); return MersenneDivisorCycles.GetCycle(modulus);
} }
#if DEBUG #if DEBUG
Console.WriteLine("Trying heuristic. Wide prime order calculation"); Console.WriteLine("Trying heuristic. Wide prime order calculation");
#endif #endif
UInt128 order = PrimeOrderCalculator.Calculate( UInt128 order = PrimeOrderCalculator.Calculate(
modulus, modulus,
_previousWidePrimeOrder, _previousWidePrimeOrder,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
if (order != UInt128.Zero) if (order != UInt128.Zero)
{ {
_previousWidePrimeOrder = order; _previousWidePrimeOrder = order;
return order; return order;
} }
_previousWidePrimeOrder = null; _previousWidePrimeOrder = null;
Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation"); Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
return MersenneDivisorCycles.GetCycle(modulus); return MersenneDivisorCycles.GetCycle(modulus);
} }
private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
private static ulong ComputeMontgomeryNPrime(ulong modulus) private static ulong ComputeMontgomeryNPrime(ulong modulus)
{ {
ulong inv = 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);
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); return unchecked(0UL - inv);
} }
} }