diff --git a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs index 090c163..238943e 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs @@ -50,7 +50,7 @@ public class GpuPow2ModBenchmarks return Pow2ModSingleBit(Input.Exponent, Input.Modulus); } - private static GpuUInt128 Pow2ModSingleBit(ulong exponent, in GpuUInt128 modulus) + private static GpuUInt128 Pow2ModSingleBit(ulong exponent, GpuUInt128 modulus) { if (modulus.IsZero || modulus == GpuUInt128.One) { diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index 48c2050..2504a3e 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -52,14 +52,31 @@ 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 Baseline() + public ulong BaselineCpu() { GetData(out var exponents, out var divisors, out _); ulong checksum = 0UL; for (int i = 0; i < SampleCount; i++) { - checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false); + checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false); + } + + return checksum; + } + + /// + /// Baseline right-to-left Montgomery ladder; measured ??.?? μs on the small sample set and ??.?? μs on the large one. + /// + [Benchmark] + public ulong BaselineGpu() + { + GetData(out var exponents, out var divisors, out _); + ulong checksum = 0UL; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false); } return checksum; @@ -69,14 +86,31 @@ 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 BaselineWithKnownCycle() + public ulong BaselineWithKnownCycleCpu() { GetData(out var exponents, out var divisors, out var cycles); ulong checksum = 0UL; for (int i = 0; i < SampleCount; i++) { - checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]); + checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]); + } + + return checksum; + } + + /// + /// 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() + { + GetData(out var exponents, out var divisors, out var cycles); + ulong checksum = 0UL; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]); } return checksum; diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs index 7bee614..f4f8aef 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs @@ -89,14 +89,32 @@ public class Pow2MontgomeryModCycleComputationBenchmarks /// large set. /// // [Benchmark(Baseline = true)] - public ulong MontgomeryWithoutCycle() + public ulong MontgomeryWithoutCycleCpu() { GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); ulong checksum = 0UL; for (int i = 0; i < SampleCount; i++) { - checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false); + checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false); + } + + return checksum; + } + + /// + /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the + /// large set. + /// + // [Benchmark(Baseline = true)] + public ulong MontgomeryWithoutCycleGpu() + { + GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); + ulong checksum = 0UL; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false); } return checksum; @@ -107,14 +125,32 @@ public class Pow2MontgomeryModCycleComputationBenchmarks /// on the large one (0.08×). /// [Benchmark] - public ulong MontgomeryWithPrecomputedCycle() + public ulong MontgomeryWithPrecomputedCycleCpu() { 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]); + checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]); + } + + return checksum; + } + + /// + /// 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×). + /// + [Benchmark] + public ulong MontgomeryWithPrecomputedCycleGpu() + { + GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles); + ulong checksum = 0UL; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]); } return checksum; @@ -125,7 +161,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline). /// // [Benchmark] - public ulong MontgomeryWithGpuCycleComputation() + public ulong MontgomeryWithGpuCycleComputationCpu() { GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); ulong checksum = 0UL; @@ -133,15 +169,33 @@ public class Pow2MontgomeryModCycleComputationBenchmarks for (int i = 0; i < SampleCount; i++) { ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus); - checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]); + checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]); } return checksum; } + /// + /// 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). + /// + // [Benchmark] + public ulong MontgomeryWithGpuCycleComputationGpu() + { + 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].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]); + } + + return checksum; + } [Benchmark] - public ulong MontgomeryWithHeuristicCycleComputation() + public ulong MontgomeryWithHeuristicCycleComputationCpu() { GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _); ulong checksum = 0UL; @@ -149,7 +203,22 @@ public class Pow2MontgomeryModCycleComputationBenchmarks for (int i = 0; i < SampleCount; i++) { ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus); - checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]); + checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]); + } + + return checksum; + } + + [Benchmark] + public ulong MontgomeryWithHeuristicCycleComputationGpu() + { + GetData(out ulong[] exponents, out MontgomeryDivisorData[] 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; diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs index 440cd96..d3ad626 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs @@ -176,7 +176,7 @@ public class Pow2MontgomeryModUInt128Benchmarks return (UInt128)result; } - private static bool ShouldUseSingleBit(in GpuUInt128 exponent) + private static bool ShouldUseSingleBit(GpuUInt128 exponent) { if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold) { @@ -257,7 +257,7 @@ public class Pow2MontgomeryModUInt128Benchmarks return result; } - private static bool IsBitSet(in GpuUInt128 value, int bitIndex) + private static bool IsBitSet(GpuUInt128 value, int bitIndex) { if (bitIndex >= 64) { @@ -267,7 +267,7 @@ public class Pow2MontgomeryModUInt128Benchmarks return ((value.Low >> bitIndex) & 1UL) != 0UL; } - private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount) + private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) { if (windowStart != 0) { diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 41da801..101f215 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,15 +1,12 @@ -using System; using System.Collections.Concurrent; -using System.Collections.Generic; using System.Numerics; -using PerfectNumbers.Core; namespace PerfectNumbers.Core.Cpu; public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester { private readonly object _sync = new(); - private readonly ConcurrentBag _sessionPool = new(); + private readonly ConcurrentBag _sessionPool = []; private ulong _divisorLimit; private ulong _lastStatusDivisor; private bool _isConfigured; @@ -213,7 +210,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong mersenne; if (prime <= 63UL) { - mersenne = (1UL << (int)prime) - 1UL; + mersenne = (1UL << (int)prime) - 1UL; hasExactMersenne = true; } else if (prime == 64UL) @@ -249,7 +246,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv else { factorCache ??= new Dictionary(8); - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) { divisorCycle = cycleCache.GetCycleLength(candidate); } @@ -287,12 +284,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv UInt128 incrementAmount = step * increments; divisor += incrementAmount; - remainder10 = AdvanceRemainder(remainder10, step10, (byte)10, increments); - remainder8 = AdvanceRemainder(remainder8, step8, (byte)8, increments); - remainder5 = AdvanceRemainder(remainder5, step5, (byte)5, increments); - remainder3 = AdvanceRemainder(remainder3, step3, (byte)3, increments); - remainder7 = AdvanceRemainder(remainder7, step7, (byte)7, increments); - remainder11 = AdvanceRemainder(remainder11, step11, (byte)11, increments); + remainder10 = AdvanceRemainder(remainder10, step10, 10, increments); + remainder8 = AdvanceRemainder(remainder8, step8, 8, increments); + remainder5 = AdvanceRemainder(remainder5, step5, 5, increments); + remainder3 = AdvanceRemainder(remainder3, step3, 3, increments); + remainder7 = AdvanceRemainder(remainder7, step7, 7, increments); + remainder11 = AdvanceRemainder(remainder11, step11, 11, increments); } processedAll = divisor > limit; @@ -301,7 +298,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { - ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData); + ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData); return residue == 1UL ? (byte)1 : (byte)0; } @@ -323,8 +320,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return remainder; } - ulong modulusValue = (ulong)modulus; - ulong scaledDelta = ((ulong)delta * (increments % modulusValue)) % modulusValue; + ulong modulusValue = modulus; + ulong scaledDelta = (delta * (increments % modulusValue)) % modulusValue; byte reducedDelta = (byte)scaledDelta; int sum = remainder + reducedDelta; if (sum >= modulus) diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs index 39fb669..4d84ec3 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepper.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs @@ -1,31 +1,20 @@ namespace PerfectNumbers.Core; -internal struct ExponentRemainderStepper +internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor) { - private readonly MontgomeryDivisorData _divisor; - private readonly ulong _modulus; - private readonly ulong _nPrime; - private readonly ulong _montgomeryOne; - private ulong _previousExponent; - private ulong _currentMontgomery; - private bool _hasState; + private readonly MontgomeryDivisorData _divisor = divisor; + private readonly ulong _modulus = divisor.Modulus; + private readonly ulong _nPrime = divisor.NPrime; + private readonly ulong _montgomeryOne = divisor.MontgomeryOne; + private ulong _previousExponent = 0UL; + private ulong _currentMontgomery = divisor.MontgomeryOne; + private bool _hasState = false; - public ExponentRemainderStepper(in MontgomeryDivisorData divisor) - { - _divisor = divisor; - _modulus = divisor.Modulus; - _nPrime = divisor.NPrime; - _montgomeryOne = divisor.MontgomeryOne; - _previousExponent = 0UL; - _currentMontgomery = divisor.MontgomeryOne; - _hasState = false; - } + public readonly bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL; - public bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL; + public readonly bool HasState => _hasState; - public bool HasState => _hasState; - - public ulong PreviousExponent => _previousExponent; + public readonly ulong PreviousExponent => _previousExponent; public void Reset() { @@ -47,7 +36,7 @@ internal struct ExponentRemainderStepper // TODO: Route these state resets through the ProcessEightBitWindows helper once the scalar // Pow2MontgomeryMod implementation adopts it so delta stepping inherits the benchmarked // 2× gains for large exponents. - _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); + _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _previousExponent = exponent; _hasState = true; return ReduceCurrent(); @@ -59,7 +48,7 @@ internal struct ExponentRemainderStepper // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks. // TODO: Replace this per-delta powmod with the upcoming windowed ladder so incremental // updates stop paying the single-bit cost highlighted in GpuPow2ModBenchmarks. - ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); + ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _previousExponent = exponent; return ReduceCurrent(); @@ -77,7 +66,7 @@ internal struct ExponentRemainderStepper { // TODO: Switch this reload to the shared windowed pow2 helper once available so CPU // residue checks align with the optimized ProcessEightBitWindows timings. - _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); + _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _previousExponent = exponent; _hasState = true; return _currentMontgomery == _montgomeryOne; @@ -88,7 +77,7 @@ internal struct ExponentRemainderStepper // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor. // TODO: Use the windowed delta pow2 helper here as well to avoid the single-bit ladder that // currently lags behind the benchmarked implementation for huge divisor cycles. - ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); + ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _previousExponent = exponent; return _currentMontgomery == _montgomeryOne; @@ -123,7 +112,7 @@ internal struct ExponentRemainderStepper { // TODO: Replace this fallback path with the upcoming ProcessEightBitWindows helper so // fresh Montgomery states also benefit from the faster pow2 ladder measured on CPUs. - _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); + _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _previousExponent = exponent; _hasState = true; isUnity = _currentMontgomery == _montgomeryOne; diff --git a/PerfectNumbers.Core/FactorEntry.cs b/PerfectNumbers.Core/FactorEntry.cs new file mode 100644 index 0000000..e4ba3be --- /dev/null +++ b/PerfectNumbers.Core/FactorEntry.cs @@ -0,0 +1,8 @@ +namespace PerfectNumbers.Core +{ + public readonly struct FactorEntry(ulong value, int exponent) + { + public readonly ulong Value = value; + public readonly int Exponent = exponent; + } +} \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index cee3beb..e3e64f8 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -267,7 +267,7 @@ public struct GpuUInt128 : IComparable, IEquatable public void Add(in ReadOnlyGpuUInt128 other) => AddInternal(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Add(in GpuUInt128 other) => AddInternal(other.High, other.Low); + public void Add(GpuUInt128 other) => AddInternal(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] private void AddInternal(ulong otherHigh, ulong otherLow) @@ -294,17 +294,17 @@ public struct GpuUInt128 : IComparable, IEquatable public void AddMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void AddMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus); + public void AddMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void AddMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus) + public void AddMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus) { ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly(); AddModInternal(value.High, value.Low, in readOnlyMod); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void AddMod(in GpuUInt128 value, in GpuUInt128 modulus) + public void AddMod(GpuUInt128 value, GpuUInt128 modulus) { ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly(); AddModInternal(value.High, value.Low, in readOnlyMod); @@ -390,7 +390,7 @@ public struct GpuUInt128 : IComparable, IEquatable } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void AddMod(in GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus); + public void AddMod(GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus); [MethodImpl(MethodImplOptions.AggressiveInlining)] public void AddMod(ulong value, ulong modulus) @@ -414,7 +414,7 @@ public struct GpuUInt128 : IComparable, IEquatable public void Sub(in ReadOnlyGpuUInt128 other) => SubInternal(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Sub(in GpuUInt128 other) => SubInternal(other.High, other.Low); + public void Sub(GpuUInt128 other) => SubInternal(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] private void SubInternal(ulong otherHigh, ulong otherLow) @@ -461,7 +461,7 @@ public struct GpuUInt128 : IComparable, IEquatable public void Mul64(in ReadOnlyGpuUInt128 other) => Mul64Internal(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Mul64(in GpuUInt128 other) => Mul64Internal(other.High, other.Low); + public void Mul64(GpuUInt128 other) => Mul64Internal(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] private void Mul64Internal(ulong otherHigh, ulong otherLow) @@ -588,7 +588,7 @@ public struct GpuUInt128 : IComparable, IEquatable } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static int TrailingZeroCount(in GpuUInt128 value) + public static int TrailingZeroCount(GpuUInt128 value) { if (value.Low != 0UL) { @@ -646,17 +646,17 @@ public struct GpuUInt128 : IComparable, IEquatable public void SubMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void SubMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus); + public void SubMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void SubMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus) + public void SubMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus) { ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly(); SubModInternal(value.High, value.Low, in readOnlyMod); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void SubMod(in GpuUInt128 value, in GpuUInt128 modulus) + public void SubMod(GpuUInt128 value, GpuUInt128 modulus) { ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly(); SubModInternal(value.High, value.Low, in readOnlyMod); @@ -727,7 +727,7 @@ public struct GpuUInt128 : IComparable, IEquatable } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Mul(in GpuUInt128 other) + public void Mul(GpuUInt128 other) { ReadOnlyGpuUInt128 readOnlyThis = AsReadOnly(); ReadOnlyGpuUInt128 readOnlyOther = other.AsReadOnly(); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 5a41f3f..2220711 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -551,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; + hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -770,7 +770,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong exponent = exponents[index]; - results[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: true); + results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index cb920cc..d38c0cc 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -811,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = phi / factor; - if (reduced.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL) + if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL) { return false; } @@ -1595,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) { - return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL; + return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL; } private static ulong CalculateByDoublingKernel(ulong prime) @@ -1795,7 +1795,7 @@ internal static partial class PrimeOrderGpuHeuristics return 0UL; } - return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false); + return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false); } private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) @@ -1811,7 +1811,7 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { ulong exponent = exponents[index]; - remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false); + remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/HeuristicFailureLog.cs b/PerfectNumbers.Core/HeuristicFailureLog.cs new file mode 100644 index 0000000..aa1b48d --- /dev/null +++ b/PerfectNumbers.Core/HeuristicFailureLog.cs @@ -0,0 +1,98 @@ +using System.Collections.Concurrent; +using System.Globalization; +using System.Text; + +namespace PerfectNumbers.Core +{ + internal static class HeuristicFailureLog + { + private const string LogFileName = "prime-order-heuristic-fallbacks.log"; + private const int FlushThreshold = 4096; + private const int MaxBufferedEntries = 16; + private static readonly ThreadLocal s_threadBuilder = new(() => StringBuilderPool.Rent(), trackAllValues: true); + private static readonly ThreadLocal s_threadEntryCount = new(() => 0, trackAllValues: true); + private static readonly object s_fileLock = new(); + private static readonly string s_logPath = Path.Combine(AppContext.BaseDirectory, LogFileName); + private static readonly UTF8Encoding s_utf8NoBom = new(false); + + static HeuristicFailureLog() + { + AppDomain.CurrentDomain.ProcessExit += static (_, _) => FlushAll(force: true); + } + + public static void Record(ulong prime, ulong? candidateOrder, HeuristicFailureReason reason) + { + string primeText = prime.ToString(CultureInfo.InvariantCulture); + string? candidateText = candidateOrder.HasValue ? candidateOrder.Value.ToString(CultureInfo.InvariantCulture) : null; + RecordInternal(primeText, candidateText, reason); + } + + public static void Record(UInt128 prime, UInt128? candidateOrder, HeuristicFailureReason reason) + { + string primeText = prime.ToString(CultureInfo.InvariantCulture); + string? candidateText = candidateOrder.HasValue ? candidateOrder.Value.ToString(CultureInfo.InvariantCulture) : null; + RecordInternal(primeText, candidateText, reason); + } + + private static void RecordInternal(string primeText, string? candidateText, HeuristicFailureReason reason) + { + StringBuilder? existing = s_threadBuilder.Value; + StringBuilder builder = existing ?? StringBuilderPool.Rent(); + int count = s_threadEntryCount.Value + 1; + bool forceFlush = count >= MaxBufferedEntries; + builder.Append(DateTime.UtcNow.ToString("o", CultureInfo.InvariantCulture)); + builder.Append(" | prime="); + builder.Append(primeText); + builder.Append(" | reason="); + builder.Append(reason); + if (candidateText is not null) + { + builder.Append(" | candidateOrder="); + builder.Append(candidateText); + } + + builder.AppendLine(); + FlushBuilderIfNeeded(builder, forceFlush); + if (builder.Length == 0) + { + count = 0; + } + + s_threadEntryCount.Value = count; + s_threadBuilder.Value = builder; + } + + private static void FlushAll(bool force) + { + foreach (StringBuilder? builder in s_threadBuilder.Values) + { + if (builder is not null) + { + FlushBuilderIfNeeded(builder, force); + } + } + } + + private static void FlushBuilderIfNeeded(StringBuilder builder, bool force) + { + if (builder.Length == 0) + { + return; + } + + if (!force && builder.Length < FlushThreshold) + { + return; + } + + string text = builder.ToString(); + builder.Clear(); + lock (s_fileLock) + { + using FileStream stream = new FileStream(s_logPath, FileMode.Append, FileAccess.Write, FileShare.ReadWrite); + using StreamWriter writer = new StreamWriter(stream, s_utf8NoBom); + writer.Write(text); + } + } + } +} \ No newline at end of file diff --git a/PerfectNumbers.Core/HeuristicFailureReason.cs b/PerfectNumbers.Core/HeuristicFailureReason.cs new file mode 100644 index 0000000..6c4c4b7 --- /dev/null +++ b/PerfectNumbers.Core/HeuristicFailureReason.cs @@ -0,0 +1,11 @@ +namespace PerfectNumbers.Core +{ + internal enum HeuristicFailureReason + { + PhiPartialFactorizationFailed, + StrictModeRequested, + HeuristicPipelineUnresolved, + HeuristicPipelineAssumedFullOrder, + } + +} \ No newline at end of file diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 1d1b1f3..d8eb5ed 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -54,15 +54,15 @@ public class MersenneDivisorCycles return snapshot; } - public static bool CycleEqualsExponent(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) + public static bool CycleEqualsExponentCpu(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) { Dictionary? cache = null; - return CycleEqualsExponent(divisor, exponent, divisorData, ref cache); + return CycleEqualsExponentCpu(divisor, exponent, divisorData, ref cache); } - public static bool CycleEqualsExponent(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary? factorCache) + public static bool CycleEqualsExponentCpu(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary? factorCache) { - if (CycleEqualsExponentForMersenneCandidate(divisor, divisorData, exponent)) + if (CycleEqualsExponentForMersenneCandidateCpu(divisor, divisorData, exponent)) { return true; } @@ -97,15 +97,15 @@ public class MersenneDivisorCycles factorCache = cache; } - if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL) + if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL) { return cycleLength == exponent; } - return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; } - public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) + public static bool CycleEqualsExponentForMersenneCandidateCpu(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) { if (!IsValidMersenneDivisorCandidate(divisor, exponent)) { @@ -126,7 +126,7 @@ public class MersenneDivisorCycles } } - return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; } private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) @@ -264,7 +264,7 @@ public class MersenneDivisorCycles return order; } - public static bool TryCalculateCycleLengthForExponent( + public static bool TryCalculateCycleLengthForExponentCpu( ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, @@ -327,7 +327,7 @@ public class MersenneDivisorCycles return false; } - cycleLength = ReduceOrder(divisorData, phi, factorCounts); + cycleLength = ReduceOrderCpu(divisorData, phi, factorCounts); return true; } @@ -471,7 +471,7 @@ public class MersenneDivisorCycles } } - private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) + private static ulong ReduceOrderCpu(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) { if (factorCounts.Count == 0) { @@ -505,7 +505,7 @@ public class MersenneDivisorCycles } ulong candidate = order / prime; - if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL) { order = candidate; continue; diff --git a/PerfectNumbers.Core/PartialFactorResult.cs b/PerfectNumbers.Core/PartialFactorResult.cs new file mode 100644 index 0000000..d47363b --- /dev/null +++ b/PerfectNumbers.Core/PartialFactorResult.cs @@ -0,0 +1,76 @@ +using System.Buffers; + +namespace PerfectNumbers.Core +{ + public sealed class PartialFactorResult + { + [ThreadStatic] + private static PartialFactorResult? s_poolHead; + + public FactorEntry[]? Factors; + private PartialFactorResult? _next; + public ulong Cofactor; + public bool FullyFactored; + public int Count; + + private PartialFactorResult() + { + } + + public static PartialFactorResult Rent(in FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count) + { + PartialFactorResult? instance = s_poolHead; + if (instance is null) + { + instance = new(); + } + else + { + s_poolHead = instance._next; + } + + instance._next = null; + instance.Factors = factors; + instance.Cofactor = cofactor; + instance.FullyFactored = fullyFactored; + instance.Count = count; + return instance; + } + + public static readonly PartialFactorResult Empty = new() { Cofactor = 1UL, FullyFactored = true }; + + public PartialFactorResult WithAdditionalPrime(ulong prime) + { + FactorEntry[]? source = Factors; + if (source is null || Count == 0) + { + FactorEntry[] local = ArrayPool.Shared.Rent(1); + local[0] = new FactorEntry(prime, 1); + return Rent(local, 1UL, true, 1); + } + + FactorEntry[] extended = ArrayPool.Shared.Rent(Count + 1); + Array.Copy(source, 0, extended, 0, Count); + extended[Count] = new FactorEntry(prime, 1); + Span span = extended.AsSpan(0, Count + 1); + span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + return Rent(extended, 1UL, true, Count + 1); + } + + public void Dispose() + { + if (this != Empty) + { + FactorEntry[]? factors = Factors; + if (factors is not null) + { + Factors = null; + ArrayPool.Shared.Return(factors, clearArray: false); + } + + _next = s_poolHead; + s_poolHead = this; + } + } + } +} \ No newline at end of file diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 4eb6e93..3e37bd3 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -8,6 +8,7 @@ public static class PerfectNumberConstants public const int MaxQForDivisorCycles = 4_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxOddPowersCount = 128; + public const int PooledArrayThreshold = 64; // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal // divisor-cycle datasets we generate offline. // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 5124839..f96168c 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -1,60 +1,31 @@ -using System; using System.Buffers; using System.Buffers.Binary; -using System.Collections.Concurrent; -using System.Collections.Generic; using System.Diagnostics; -using System.Globalization; -using System.IO; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Security.Cryptography; -using System.Text; using PerfectNumbers.Core.Gpu; -using System.Threading; namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - - internal enum PrimeOrderMode { Heuristic, Strict, } - private enum HeuristicFailureReason - { - PhiPartialFactorizationFailed, - StrictModeRequested, - HeuristicPipelineUnresolved, - HeuristicPipelineAssumedFullOrder, - } - - internal readonly struct PrimeOrderSearchConfig - { - public PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) - { - SmallFactorLimit = smallFactorLimit; - PollardRhoMilliseconds = pollardRhoMilliseconds; - MaxPowChecks = maxPowChecks; - Mode = mode; - } - - public uint SmallFactorLimit { get; } - - public int PollardRhoMilliseconds { get; } - - public int MaxPowChecks { get; } - - public PrimeOrderMode Mode { get; } - - public static PrimeOrderSearchConfig HeuristicDefault => new(smallFactorLimit: 100_000, pollardRhoMilliseconds: 128, maxPowChecks: 24, PrimeOrderMode.Heuristic); - + internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) + { + public static PrimeOrderSearchConfig HeuristicDefault => new(smallFactorLimit: 100_000, pollardRhoMilliseconds: 128, maxPowChecks: 24, PrimeOrderMode.Heuristic); public static PrimeOrderSearchConfig StrictDefault => new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); + + public readonly uint SmallFactorLimit = smallFactorLimit; + public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; + public readonly int MaxPowChecks = maxPowChecks; + public readonly PrimeOrderMode Mode = mode; } internal enum PrimeOrderHeuristicDevice @@ -88,29 +59,30 @@ internal static partial class PrimeOrderCalculator private static Dictionary? s_partialFactorCounts; private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; - private readonly struct Pow2ModeScope - { - private readonly bool _previousInitialized; - private readonly bool _previousAllow; - private readonly PrimeOrderHeuristicDevice _previousDevice; - public Pow2ModeScope(PrimeOrderHeuristicDevice device) - { - _previousInitialized = s_pow2ModeInitialized; - _previousAllow = s_allowGpuPow2; - _previousDevice = s_deviceMode; - s_pow2ModeInitialized = true; - s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; - s_deviceMode = device; - } + private readonly struct Pow2ModeScope + { + private readonly bool _previousInitialized; + private readonly bool _previousAllow; + private readonly PrimeOrderHeuristicDevice _previousDevice; - public void Dispose() - { - s_allowGpuPow2 = _previousAllow; - s_pow2ModeInitialized = _previousInitialized; - s_deviceMode = _previousDevice; - } - } + public Pow2ModeScope(PrimeOrderHeuristicDevice device) + { + _previousInitialized = s_pow2ModeInitialized; + _previousAllow = s_allowGpuPow2; + _previousDevice = s_deviceMode; + s_pow2ModeInitialized = true; + s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; + s_deviceMode = device; + } + + public void Dispose() + { + s_allowGpuPow2 = _previousAllow; + s_pow2ModeInitialized = _previousInitialized; + s_deviceMode = _previousDevice; + } + } private readonly struct DebugLoggingScope { @@ -184,17 +156,22 @@ internal static partial class PrimeOrderCalculator return gpuOrder; } - using PartialFactorResult phiFactors = PartialFactor(phi, config); + PartialFactorResult phiFactors = PartialFactor(phi, config); - if (phiFactors.Factors is null) - { - return FinishStrictly(prime, divisorData); - } + ulong result; + if (phiFactors.Factors is null) + { + result = FinishStrictlyCpu(prime, divisorData); + phiFactors.Dispose(); + return result; + } - return RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors); + result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); + phiFactors.Dispose(); + return result; } - private static ulong RunHeuristicPipeline( + private static ulong RunHeuristicPipelineCpu( ulong prime, ulong? previousOrder, in PrimeOrderSearchConfig config, @@ -202,25 +179,25 @@ internal static partial class PrimeOrderCalculator ulong phi, PartialFactorResult phiFactors) { - if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData)) + if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) { return phi; } - ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData); - candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData); + ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); + candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - if (TryConfirmOrder(prime, candidateOrder, divisorData, config)) + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config)) { return candidateOrder; } if (config.Mode == PrimeOrderMode.Strict) { - return FinishStrictly(prime, divisorData); + return FinishStrictlyCpu(prime, divisorData); } - if (TryHeuristicFinish(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) + if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) { return order; } @@ -228,12 +205,12 @@ internal static partial class PrimeOrderCalculator return candidateOrder; } - private static ulong FinishStrictly(ulong prime, in MontgomeryDivisorData divisorData) + private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData) { - return CalculateByFactorization(prime, divisorData); + return CalculateByFactorizationCpu(prime, divisorData); } - private static bool TrySpecialMax(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) { ReadOnlySpan factorSpan = factors.Factors; int length = factors.Count; @@ -241,7 +218,7 @@ internal static partial class PrimeOrderCalculator { ulong factor = factorSpan[i].Value; ulong reduced = phi / factor; - if (Pow2EqualsOne(reduced, prime, divisorData)) + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) { return false; } @@ -250,13 +227,13 @@ internal static partial class PrimeOrderCalculator return true; } - private static ulong InitializeStartingOrder(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) + private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) { ulong order = phi; if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) { ulong half = phi >> 1; - if (Pow2EqualsOne(half, prime, divisorData)) + if (Pow2EqualsOneCpu(half, prime, divisorData)) { order = half; } @@ -265,7 +242,7 @@ internal static partial class PrimeOrderCalculator return order; } - private static ulong ExponentLowering(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) { FactorEntry[]? tempArray = null; try @@ -298,7 +275,7 @@ internal static partial class PrimeOrderCalculator } ulong reduced = order / primeFactor; - if (Pow2EqualsOne(reduced, prime, divisorData)) + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) { order = reduced; continue; @@ -319,7 +296,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmOrder(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { if (order == 0UL) { @@ -328,7 +305,7 @@ internal static partial class PrimeOrderCalculator // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`. // DebugLog("Verifying a^order ≡ 1 (mod p)"); - if (!Pow2EqualsOne(order, prime, divisorData)) + if (!Pow2EqualsOneCpu(order, prime, divisorData)) { return false; } @@ -380,7 +357,7 @@ internal static partial class PrimeOrderCalculator } reduced /= primeFactor; - if (Pow2EqualsOne(reduced, prime, divisorData)) + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) { return false; } @@ -395,7 +372,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryHeuristicFinish( + private static bool TryHeuristicFinishCpu( ulong prime, ulong order, ulong? previousOrder, @@ -528,14 +505,14 @@ internal static partial class PrimeOrderCalculator } else { - equalsOne = Pow2EqualsOne(candidate, prime, divisorData); + equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData); } if (!equalsOne) { continue; } - if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget)) + if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) { continue; } @@ -733,7 +710,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmCandidate(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) + private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) { PartialFactorResult factorization = PartialFactor(candidate, config); try @@ -782,7 +759,7 @@ internal static partial class PrimeOrderCalculator } powUsed++; - if (Pow2EqualsOne(reduced, prime, divisorData)) + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) { return false; } @@ -797,7 +774,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) + private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) { if (IsGpuPow2Allowed) { @@ -808,111 +785,7 @@ internal static partial class PrimeOrderCalculator } } - return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; - } - - private static class HeuristicFailureLog - { - private const string LogFileName = "prime-order-heuristic-fallbacks.log"; - private const int FlushThreshold = 4096; - private const int MaxBufferedEntries = 16; - private static readonly ConcurrentQueue s_builderPool = new(); - private static readonly ThreadLocal s_threadBuilder = new(() => AcquireBuilder(), trackAllValues: true); - private static readonly ThreadLocal s_threadEntryCount = new(() => 0, trackAllValues: true); - private static readonly object s_fileLock = new(); - private static readonly string s_logPath = Path.Combine(AppContext.BaseDirectory, LogFileName); - private static readonly UTF8Encoding s_utf8NoBom = new(false); - - static HeuristicFailureLog() - { - AppDomain.CurrentDomain.ProcessExit += static (_, _) => FlushAll(force: true); - } - - public static void Record(ulong prime, ulong? candidateOrder, HeuristicFailureReason reason) - { - string primeText = prime.ToString(CultureInfo.InvariantCulture); - string? candidateText = candidateOrder.HasValue ? candidateOrder.Value.ToString(CultureInfo.InvariantCulture) : null; - RecordInternal(primeText, candidateText, reason); - } - - public static void Record(UInt128 prime, UInt128? candidateOrder, HeuristicFailureReason reason) - { - string primeText = prime.ToString(CultureInfo.InvariantCulture); - string? candidateText = candidateOrder.HasValue ? candidateOrder.Value.ToString(CultureInfo.InvariantCulture) : null; - RecordInternal(primeText, candidateText, reason); - } - - private static void RecordInternal(string primeText, string? candidateText, HeuristicFailureReason reason) - { - StringBuilder? existing = s_threadBuilder.Value; - StringBuilder builder = existing ?? AcquireBuilder(); - int count = s_threadEntryCount.Value + 1; - bool forceFlush = count >= MaxBufferedEntries; - builder.Append(DateTime.UtcNow.ToString("o", CultureInfo.InvariantCulture)); - builder.Append(" | prime="); - builder.Append(primeText); - builder.Append(" | reason="); - builder.Append(reason); - if (candidateText is not null) - { - builder.Append(" | candidateOrder="); - builder.Append(candidateText); - } - - builder.AppendLine(); - FlushBuilderIfNeeded(builder, forceFlush); - if (builder.Length == 0) - { - count = 0; - } - - s_threadEntryCount.Value = count; - s_threadBuilder.Value = builder; - } - - private static void FlushAll(bool force) - { - foreach (StringBuilder? builder in s_threadBuilder.Values) - { - if (builder is not null) - { - FlushBuilderIfNeeded(builder, force); - } - } - } - - private static void FlushBuilderIfNeeded(StringBuilder builder, bool force) - { - if (builder.Length == 0) - { - return; - } - - if (!force && builder.Length < FlushThreshold) - { - return; - } - - string text = builder.ToString(); - builder.Clear(); - lock (s_fileLock) - { - using FileStream stream = new FileStream(s_logPath, FileMode.Append, FileAccess.Write, FileShare.ReadWrite); - using StreamWriter writer = new StreamWriter(stream, s_utf8NoBom); - writer.Write(text); - } - } - - private static StringBuilder AcquireBuilder() - { - if (s_builderPool.TryDequeue(out StringBuilder? builder)) - { - builder.Clear(); - return builder; - } - - return new StringBuilder(FlushThreshold); - } + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; } private static List AcquirePendingCompositeList(int capacityHint) @@ -952,16 +825,15 @@ internal static partial class PrimeOrderCalculator private static Dictionary AcquireFactorCountDictionary(int capacityHint) { Dictionary? dictionary = s_partialFactorCounts; - int adjustedCapacity = Math.Max(capacityHint, 8); if (dictionary is null) { - dictionary = new Dictionary(adjustedCapacity); + dictionary = new Dictionary(capacityHint); s_partialFactorCounts = dictionary; return dictionary; } dictionary.Clear(); - dictionary.EnsureCapacity(adjustedCapacity); + dictionary.EnsureCapacity(capacityHint); return dictionary; } @@ -969,7 +841,7 @@ internal static partial class PrimeOrderCalculator { if (value <= 1UL) { - return PartialFactorResult.RentEmpty(); + return PartialFactorResult.Empty; } const int FactorSlotCount = GpuSmallPrimeFactorSlots; @@ -1459,7 +1331,7 @@ internal static partial class PrimeOrderCalculator } } - private static ulong CalculateByFactorization(ulong prime, in MontgomeryDivisorData divisorData) + private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) { ulong phi = prime - 1UL; Dictionary counts = new(capacity: 8); @@ -1486,7 +1358,7 @@ internal static partial class PrimeOrderCalculator } ulong candidate = order / primeFactor; - if (Pow2EqualsOne(candidate, prime, divisorData)) + if (Pow2EqualsOneCpu(candidate, prime, divisorData)) { order = candidate; continue; @@ -1557,96 +1429,4 @@ internal static partial class PrimeOrderCalculator public readonly long Secondary = secondary; public readonly long Tertiary = tertiary; } - - private readonly struct FactorEntry(ulong value, int exponent) - { - public readonly ulong Value = value; - public readonly int Exponent = exponent; - } - - private sealed class PartialFactorResult : IDisposable - { - [ThreadStatic] - private static PartialFactorResult? s_poolHead; - - private PartialFactorResult? _next; - private bool _returned; - - public FactorEntry[]? Factors; - public ulong Cofactor; - public bool FullyFactored; - public int Count; - - private PartialFactorResult() - { - } - - public static PartialFactorResult Rent(FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count) - { - PartialFactorResult? instance = s_poolHead; - if (instance is null) - { - instance = new PartialFactorResult(); - } - else - { - s_poolHead = instance._next; - } - - instance._next = null; - instance._returned = false; - instance.Factors = factors; - instance.Cofactor = cofactor; - instance.FullyFactored = fullyFactored; - instance.Count = count; - return instance; - } - - public static PartialFactorResult RentEmpty() - { - return Rent(null, 1UL, true, 0); - } - - public PartialFactorResult WithAdditionalPrime(ulong prime) - { - FactorEntry[]? source = Factors; - if (source is null || Count == 0) - { - FactorEntry[] local = ArrayPool.Shared.Rent(1); - local[0] = new FactorEntry(prime, 1); - return Rent(local, 1UL, true, 1); - } - - FactorEntry[] extended = ArrayPool.Shared.Rent(Count + 1); - Array.Copy(source, 0, extended, 0, Count); - extended[Count] = new FactorEntry(prime, 1); - Span span = extended.AsSpan(0, Count + 1); - span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - return Rent(extended, 1UL, true, Count + 1); - } - - public void Dispose() - { - if (_returned) - { - return; - } - - _returned = true; - - FactorEntry[]? factors = Factors; - if (factors is not null) - { - Factors = null; - ArrayPool.Shared.Return(factors, clearArray: false); - } - - Cofactor = 0UL; - FullyFactored = false; - Count = 0; - - _next = s_poolHead; - s_poolHead = this; - } - } } diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 7f62518..e7bf9c3 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -218,7 +218,7 @@ public static class UInt128Extensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool IsBitSet(in GpuUInt128 value, int bitIndex) + private static bool IsBitSet(GpuUInt128 value, int bitIndex) { if (bitIndex >= 64) { @@ -229,18 +229,16 @@ public static class UInt128Extensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount) + private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) { + ulong mask = (1UL << windowBitCount) - 1UL; if (windowStart != 0) { - GpuUInt128 shifted = exponent; - shifted.ShiftRight(windowStart); - ulong mask = (1UL << windowBitCount) - 1UL; - return shifted.Low & mask; + exponent.ShiftRight(windowStart); + return exponent.Low & mask; } - ulong directMask = (1UL << windowBitCount) - 1UL; - return exponent.Low & directMask; + return exponent.Low & mask; } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index f46b61c..f9ae011 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -1,4 +1,5 @@ using System; +using System.Buffers; using System.Numerics; using System.Runtime.CompilerServices; using ILGPU.Algorithms; @@ -76,7 +77,7 @@ public static class ULongExtensions public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift) { ulong high = value >> shift; - ulong mask = (1UL << shift) - 1UL; + ulong mask = (1UL << shift) - 1UL; ulong low = value & mask; UInt128 highContribution = (UInt128)high * multiplier; @@ -326,7 +327,7 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { ulong modulus = divisor.Modulus; if (exponent == 0UL) @@ -345,7 +346,78 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount); + ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount); + + int index = bitLength - 1; + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MontgomeryMultiply(multiplier, modulus, nPrime); + + index = windowStart - 1; + } + + if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold) + { + ArrayPool.Shared.Return(oddPowers); + } + + if (keepMontgomery) + { + return result; + } + + return result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + ulong modulus = divisor.Modulus; + if (exponent == 0UL) + { + return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + } + + if (exponent <= Pow2WindowFallbackThreshold) + { + return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + ulong[] oddPowers = InitializeMontgomeryOddPowersGpu(divisor, modulus, nPrime, oddPowerCount); int index = bitLength - 1; while (index >= 0) @@ -461,7 +533,30 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) + private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) + { + ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold ? new ulong[oddPowerCount] : ArrayPool.Shared.Rent(oddPowerCount); + + oddPowers[0] = divisor.MontgomeryTwo; + if (oddPowerCount == 1) + { + return oddPowers; + } + + ulong previous, + square = divisor.MontgomeryTwoSquared; + + for (int i = 1; i < oddPowerCount; i++) + { + previous = oddPowers[i - 1]; + oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); + } + + return oddPowers; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) { ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount]; oddPowers[0] = divisor.MontgomeryTwo; @@ -481,17 +576,31 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) + public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) { ulong rotationCount = exponent % cycleLength; - return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false); + return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor) + public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) { - return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false); + ulong rotationCount = exponent % cycleLength; + return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false); } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus) {