From c6ad53fffbd2db59c656b517711c99f46bbbac11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 18 Oct 2025 12:22:08 +0200 Subject: [PATCH] Restore Pow2 benchmark data reuse --- ...yModBenchmarks.cs => Pow2ModBenchmarks.cs} | 270 ++++++++++++++- .../Pow2MontgomeryModUInt128Benchmarks.cs | 308 ------------------ EvenPerfectBitScanner.Benchmarks/Program.cs | 4 +- 3 files changed, 260 insertions(+), 322 deletions(-) rename EvenPerfectBitScanner.Benchmarks/{Pow2MontgomeryModBenchmarks.cs => Pow2ModBenchmarks.cs} (80%) delete mode 100644 EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs similarity index 80% rename from EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs rename to EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs index f031f83..9584a08 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs @@ -4,6 +4,7 @@ using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using System.Numerics; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; namespace EvenPerfectBitScanner.Benchmarks; @@ -20,16 +21,18 @@ namespace EvenPerfectBitScanner.Benchmarks; /// [MemoryDiagnoser] [SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] -public class Pow2MontgomeryModBenchmarks +public class Pow2ModBenchmarks { private const int SampleCount = 256; private const ulong OneShiftedLeft63 = 1UL << 63; private const long OneShiftedLeft60 = 1L << 60; + private const ulong Pow2WindowFallbackThreshold = 32UL; private MontgomeryDivisorData[] _smallDivisors = Array.Empty(); private MontgomeryDivisorData[] _largeDivisors = Array.Empty(); private MontgomeryDivisorData[] _veryLargeDivisors = Array.Empty(); private UInt128[] _wideModuli = Array.Empty(); + private UInt128[] _wideExponents = Array.Empty(); private ulong[] _smallExponents = Array.Empty(); private ulong[] _largeExponents = Array.Empty(); private ulong[] _veryLargeExponents = Array.Empty(); @@ -58,6 +61,7 @@ public class Pow2MontgomeryModBenchmarks _largeDivisors = cache.LargeDivisors; _veryLargeDivisors = cache.VeryLargeDivisors; _wideModuli = cache.WideModuli; + _wideExponents = cache.WideExponents; _smallExponents = cache.SmallExponents; _largeExponents = cache.LargeExponents; _veryLargeExponents = cache.VeryLargeExponents; @@ -70,8 +74,8 @@ public class Pow2MontgomeryModBenchmarks /// /// Baseline right-to-left Montgomery ladder; measured 36.36 μs on the small sample set and 117.32 μs on the large one. /// - [Benchmark(Baseline = true)] - public ulong BaselineCpu() + [Benchmark] + public ulong MontgomeryWindowedCpu() { GetData(out var exponents, out var divisors, out _); ulong checksum = 0UL; @@ -88,7 +92,7 @@ public class Pow2MontgomeryModBenchmarks /// Baseline right-to-left Montgomery ladder; measured ??.?? μs on the small sample set and ??.?? μs on the large one. /// [Benchmark] - public ulong BaselineGpu() + public ulong MontgomeryWindowedGpu() { GetData(out var exponents, out var divisors, out _); ulong checksum = 0UL; @@ -105,7 +109,7 @@ public class Pow2MontgomeryModBenchmarks /// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline). /// [Benchmark] - public ulong BaselineWithKnownCycleCpu() + public ulong MontgomeryWindowedWithCycleCpu() { GetData(out var exponents, out var divisors, out var cycles); ulong checksum = 0UL; @@ -122,7 +126,7 @@ public class Pow2MontgomeryModBenchmarks /// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline). /// [Benchmark] - public ulong BaselineWithKnownCycleGpu() + public ulong MontgomeryWindowedWithCycleGpu() { GetData(out var exponents, out var divisors, out var cycles); ulong checksum = 0UL; @@ -271,8 +275,8 @@ public class Pow2MontgomeryModBenchmarks /// /// Classic binary ladder without Montgomery reduction; runs 24.75 μs on small inputs (0.68× baseline) but collapses to 295.29 μs on large ones (2.52×). /// - [Benchmark] - public ulong BinaryMod() + [Benchmark(Baseline = true)] + public ulong StandardPow2Mod() { GetData(out var exponents, out var divisors, out _); ulong checksum = 0UL; @@ -319,6 +323,38 @@ public class Pow2MontgomeryModBenchmarks return checksum; } + /// + /// Current UInt128-native windowed Montgomery ladder that avoids GpuUInt128 on CPU paths. + /// + [Benchmark] + public UInt128 MontgomeryWindowedUInt128() + { + UInt128 checksum = UInt128.Zero; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= _wideExponents[i].Pow2MontgomeryModWindowed(_wideModuli[i]); + } + + return checksum; + } + + /// + /// Previous CPU shim that reused the GPU-oriented GpuUInt128 implementation. + /// + [Benchmark] + public UInt128 MontgomeryWindowedGpuShimUInt128() + { + UInt128 checksum = UInt128.Zero; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= Pow2MontgomeryModWindowedGpuShim(_wideExponents[i], _wideModuli[i]); + } + + return checksum; + } + private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles) { switch (Scale) @@ -350,6 +386,7 @@ public class Pow2MontgomeryModBenchmarks MontgomeryDivisorData[] largeDivisors, MontgomeryDivisorData[] veryLargeDivisors, UInt128[] wideModuli, + UInt128[] wideExponents, ulong[] smallExponents, ulong[] largeExponents, ulong[] veryLargeExponents, @@ -362,6 +399,7 @@ public class Pow2MontgomeryModBenchmarks LargeDivisors = largeDivisors; VeryLargeDivisors = veryLargeDivisors; WideModuli = wideModuli; + WideExponents = wideExponents; SmallExponents = smallExponents; LargeExponents = largeExponents; VeryLargeExponents = veryLargeExponents; @@ -375,6 +413,7 @@ public class Pow2MontgomeryModBenchmarks public MontgomeryDivisorData[] LargeDivisors { get; } public MontgomeryDivisorData[] VeryLargeDivisors { get; } public UInt128[] WideModuli { get; } + public UInt128[] WideExponents { get; } public ulong[] SmallExponents { get; } public ulong[] LargeExponents { get; } public ulong[] VeryLargeExponents { get; } @@ -394,11 +433,13 @@ public class Pow2MontgomeryModBenchmarks { Random random = new(13); Random heuristicRandom = new(113); + Random wideRandom = new(19); var smallDivisors = new MontgomeryDivisorData[SampleCount]; var largeDivisors = new MontgomeryDivisorData[SampleCount]; var veryLargeDivisors = new MontgomeryDivisorData[SampleCount]; var wideModuli = new UInt128[SampleCount]; + var wideExponents = new UInt128[SampleCount]; var smallExponents = new ulong[SampleCount]; var largeExponents = new ulong[SampleCount]; var veryLargeExponents = new ulong[SampleCount]; @@ -427,8 +468,9 @@ public class Pow2MontgomeryModBenchmarks veryLargeExponents[i] = NextVeryLargeExponent(heuristicRandom); veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus, ref previousPrimeOrder); - UInt128 wideModulus = NextWideOddModulus(heuristicRandom); + UInt128 wideModulus = NextWideOddModulus(wideRandom); wideModuli[i] = wideModulus; + wideExponents[i] = NextWideExponent(wideRandom); wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus, ref previousWidePrimeOrder); } @@ -437,6 +479,7 @@ public class Pow2MontgomeryModBenchmarks largeDivisors, veryLargeDivisors, wideModuli, + wideExponents, smallExponents, largeExponents, veryLargeExponents, @@ -483,10 +526,37 @@ public class Pow2MontgomeryModBenchmarks { Span buffer = stackalloc byte[16]; random.NextBytes(buffer); + ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); - ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]) | (1UL << 63); - UInt128 value = ((UInt128)high << 64) | low; - return value | UInt128.One; + ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]); + if (high == 0UL) + { + high = 1UL; + } + + UInt128 modulus = ((UInt128)high << 64) | low; + modulus |= UInt128.One; + + if (modulus <= UInt128.One) + { + modulus += 2UL; + } + + return modulus; + } + + private static UInt128 NextWideExponent(Random random) + { + Span buffer = stackalloc byte[16]; + random.NextBytes(buffer); + ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); + ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]); + if (high == 0UL) + { + high = 1UL; + } + + return ((UInt128)high << 64) | low; } private static ulong CalculateCycleLengthWithHeuristics(ulong modulus, ref ulong? previousPrimeOrder) @@ -851,6 +921,182 @@ public class Pow2MontgomeryModBenchmarks return (ulong)result; } + private static UInt128 Pow2MontgomeryModWindowedGpuShim(UInt128 exponent, UInt128 modulus) + { + if (modulus == UInt128.One) + { + return UInt128.Zero; + } + + if (exponent == UInt128.Zero) + { + return UInt128.One % modulus; + } + + GpuUInt128 modulusGpu = new(modulus); + GpuUInt128 exponentGpu = new(exponent); + GpuUInt128 baseValue = new(2UL); + if (baseValue.CompareTo(modulusGpu) >= 0) + { + baseValue.Sub(modulusGpu); + } + + if (ShouldUseSingleBit(exponentGpu)) + { + GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponentGpu, modulusGpu, baseValue); + return (UInt128)singleBitResult; + } + + int bitLength = exponentGpu.GetBitLength(); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + Span oddPowers = stackalloc GpuUInt128[oddPowerCount]; + InitializeOddPowers(baseValue, modulusGpu, oddPowers); + + GpuUInt128 result = GpuUInt128.One; + int index = bitLength - 1; + + while (index >= 0) + { + if (!IsBitSet(exponentGpu, index)) + { + result.MulMod(result, modulusGpu); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (!IsBitSet(exponentGpu, windowStart)) + { + windowStart++; + } + + int windowBitCount = index - windowStart + 1; + for (int square = 0; square < windowBitCount; square++) + { + result.MulMod(result, modulusGpu); + } + + ulong windowValue = ExtractWindowValue(exponentGpu, windowStart, windowBitCount); + int tableIndex = (int)((windowValue - 1UL) >> 1); + GpuUInt128 factor = oddPowers[tableIndex]; + result.MulMod(factor, modulusGpu); + + index = windowStart - 1; + } + + return (UInt128)result; + } + + private static bool ShouldUseSingleBit(GpuUInt128 exponent) + { + if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold) + { + return true; + } + + return false; + } + + private static int GetWindowSize(int bitLength) + { + if (bitLength <= 8) + { + return Math.Max(bitLength, 1); + } + + if (bitLength <= 23) + { + return 4; + } + + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + return 8; + } + + private static void InitializeOddPowers(GpuUInt128 baseValue, GpuUInt128 modulus, Span oddPowers) + { + oddPowers[0] = baseValue; + if (oddPowers.Length == 1) + { + return; + } + + baseValue.MulMod(baseValue, modulus); + for (int i = 1; i < oddPowers.Length; i++) + { + oddPowers[i] = oddPowers[i - 1]; + oddPowers[i].MulMod(baseValue, modulus); + } + } + + private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) + { + GpuUInt128 result = GpuUInt128.One; + + while (!exponent.IsZero) + { + if ((exponent.Low & 1UL) != 0UL) + { + result.MulMod(baseValue, modulus); + } + + exponent.ShiftRight(1); + if (exponent.IsZero) + { + break; + } + + baseValue.MulMod(baseValue, modulus); + } + + return result; + } + + private static bool IsBitSet(GpuUInt128 value, int bitIndex) + { + if (bitIndex >= 64) + { + return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; + } + + return ((value.Low >> bitIndex) & 1UL) != 0UL; + } + + private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) + { + if (windowStart != 0) + { + GpuUInt128 shifted = exponent; + shifted.ShiftRight(windowStart); + ulong mask = (1UL << windowBitCount) - 1UL; + return shifted.Low & mask; + } + + ulong directMask = (1UL << windowBitCount) - 1UL; + return exponent.Low & directMask; + } + private static ulong MultiplyMod(ulong x, ulong y, ulong modulus) { UInt128 product = (UInt128)x * y; diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs deleted file mode 100644 index bffcaf8..0000000 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs +++ /dev/null @@ -1,308 +0,0 @@ -using System.Buffers.Binary; -using BenchmarkDotNet.Attributes; -using BenchmarkDotNet.Jobs; -using PerfectNumbers.Core; -using PerfectNumbers.Core.Gpu; - -namespace EvenPerfectBitScanner.Benchmarks; - -[MemoryDiagnoser] -[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] -public class Pow2MontgomeryModUInt128Benchmarks -{ - private const int SampleCount = 256; - private const ulong Pow2WindowFallbackThreshold = 32UL; - - private UInt128[] _moduli = Array.Empty(); - private UInt128[] _exponents = Array.Empty(); - - [GlobalSetup] - public void Setup() - { - WideSampleData cache = WideSampleDataCache.Instance; - _moduli = cache.Moduli; - _exponents = cache.Exponents; - } - - private sealed class WideSampleData - { - public WideSampleData(UInt128[] moduli, UInt128[] exponents) - { - Moduli = moduli; - Exponents = exponents; - } - - public UInt128[] Moduli { get; } - public UInt128[] Exponents { get; } - } - - private static class WideSampleDataCache - { - private static readonly Lazy Cache = new(Create); - - public static WideSampleData Instance => Cache.Value; - - private static WideSampleData Create() - { - Random random = new(19); - var moduli = new UInt128[SampleCount]; - var exponents = new UInt128[SampleCount]; - - for (int i = 0; i < SampleCount; i++) - { - moduli[i] = NextWideOddModulus(random); - exponents[i] = NextWideExponent(random); - } - - return new WideSampleData(moduli, exponents); - } - - private static UInt128 NextWideOddModulus(Random random) - { - Span buffer = stackalloc byte[16]; - random.NextBytes(buffer); - - ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); - ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]); - if (high == 0UL) - { - high = 1UL; - } - - UInt128 modulus = ((UInt128)high << 64) | low; - modulus |= UInt128.One; - - if (modulus <= UInt128.One) - { - modulus += 2UL; - } - - return modulus; - } - - private static UInt128 NextWideExponent(Random random) - { - Span buffer = stackalloc byte[16]; - random.NextBytes(buffer); - - ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); - ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]); - if (high == 0UL) - { - high = 1UL; - } - - return ((UInt128)high << 64) | low; - } - } - - /// - /// Current UInt128-native windowed Montgomery ladder that avoids GpuUInt128 on CPU paths. - /// - [Benchmark(Baseline = true)] - public UInt128 UInt128Windowed() - { - UInt128 checksum = UInt128.Zero; - - for (int i = 0; i < SampleCount; i++) - { - checksum ^= _exponents[i].Pow2MontgomeryModWindowed(_moduli[i]); - } - - return checksum; - } - - /// - /// Previous CPU shim that reused the GPU-oriented GpuUInt128 implementation. - /// - [Benchmark] - public UInt128 GpuStructShim() - { - UInt128 checksum = UInt128.Zero; - - for (int i = 0; i < SampleCount; i++) - { - checksum ^= Pow2MontgomeryModWindowedGpuShim(_exponents[i], _moduli[i]); - } - - return checksum; - } - - private static UInt128 Pow2MontgomeryModWindowedGpuShim(UInt128 exponent, UInt128 modulus) - { - if (modulus == UInt128.One) - { - return UInt128.Zero; - } - - if (exponent == UInt128.Zero) - { - return UInt128.One % modulus; - } - - GpuUInt128 modulusGpu = new(modulus); - GpuUInt128 exponentGpu = new(exponent); - GpuUInt128 baseValue = new(2UL); - if (baseValue.CompareTo(modulusGpu) >= 0) - { - baseValue.Sub(modulusGpu); - } - - if (ShouldUseSingleBit(exponentGpu)) - { - GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponentGpu, modulusGpu, baseValue); - return (UInt128)singleBitResult; - } - - int bitLength = exponentGpu.GetBitLength(); - int windowSize = GetWindowSize(bitLength); - int oddPowerCount = 1 << (windowSize - 1); - - Span oddPowers = stackalloc GpuUInt128[oddPowerCount]; - InitializeOddPowers(baseValue, modulusGpu, oddPowers); - - GpuUInt128 result = GpuUInt128.One; - int index = bitLength - 1; - - while (index >= 0) - { - if (!IsBitSet(exponentGpu, index)) - { - result.MulMod(result, modulusGpu); - index--; - continue; - } - - int windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } - - while (!IsBitSet(exponentGpu, windowStart)) - { - windowStart++; - } - - int windowBitCount = index - windowStart + 1; - for (int square = 0; square < windowBitCount; square++) - { - result.MulMod(result, modulusGpu); - } - - ulong windowValue = ExtractWindowValue(exponentGpu, windowStart, windowBitCount); - int tableIndex = (int)((windowValue - 1UL) >> 1); - GpuUInt128 factor = oddPowers[tableIndex]; - result.MulMod(factor, modulusGpu); - - index = windowStart - 1; - } - - return (UInt128)result; - } - - private static bool ShouldUseSingleBit(GpuUInt128 exponent) - { - if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold) - { - return true; - } - - return false; - } - - private static int GetWindowSize(int bitLength) - { - if (bitLength <= 8) - { - return Math.Max(bitLength, 1); - } - - if (bitLength <= 23) - { - return 4; - } - - if (bitLength <= 79) - { - return 5; - } - - if (bitLength <= 239) - { - return 6; - } - - if (bitLength <= 671) - { - return 7; - } - - return 8; - } - - private static void InitializeOddPowers(GpuUInt128 baseValue, GpuUInt128 modulus, Span oddPowers) - { - oddPowers[0] = baseValue; - if (oddPowers.Length == 1) - { - return; - } - - // Reusing baseValue to hold base^2 for the shared odd-power ladder. - baseValue.MulMod(baseValue, modulus); - for (int i = 1; i < oddPowers.Length; i++) - { - oddPowers[i] = oddPowers[i - 1]; - oddPowers[i].MulMod(baseValue, modulus); - } - } - - private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) - { - GpuUInt128 result = GpuUInt128.One; - - while (!exponent.IsZero) - { - if ((exponent.Low & 1UL) != 0UL) - { - result.MulMod(baseValue, modulus); - } - - exponent.ShiftRight(1); - if (exponent.IsZero) - { - break; - } - - // Reusing baseValue to store the squared base for the next iteration. - baseValue.MulMod(baseValue, modulus); - } - - return result; - } - - private static bool IsBitSet(GpuUInt128 value, int bitIndex) - { - if (bitIndex >= 64) - { - return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; - } - - return ((value.Low >> bitIndex) & 1UL) != 0UL; - } - - private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) - { - if (windowStart != 0) - { - GpuUInt128 shifted = exponent; - shifted.ShiftRight(windowStart); - ulong mask = (1UL << windowBitCount) - 1UL; - return shifted.Low & mask; - } - - ulong directMask = (1UL << windowBitCount) - 1UL; - return exponent.Low & directMask; - } -} diff --git a/EvenPerfectBitScanner.Benchmarks/Program.cs b/EvenPerfectBitScanner.Benchmarks/Program.cs index c42d0dd..5bdab2a 100644 --- a/EvenPerfectBitScanner.Benchmarks/Program.cs +++ b/EvenPerfectBitScanner.Benchmarks/Program.cs @@ -9,9 +9,9 @@ public static class Program // var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator); // BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config); - Pow2MontgomeryModBenchmarks benchmarks = new() + Pow2ModBenchmarks benchmarks = new() { - Scale = Pow2MontgomeryModBenchmarks.InputScale.Large + Scale = Pow2ModBenchmarks.InputScale.Large }; GpuPrimeWorkLimiter.SetLimit(1);