mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Various optimizations and cleanups. "in" removed from GpuUInt128 to avoid defensive copies from IL.

Browse Source
This commit is contained in:
mlemanczyk committed 2025-10-10 09:26:25 +02:00
1 parent 70eca5126b
commit ca0cadc8a3
18 files changed
+565 -395
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs
+1 -1
View File
@@ -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)
{
EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs
+38 -4
View File
@@ -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.
/// </summary>
[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;
}
/// <summary>
/// Baseline right-to-left Montgomery ladder; measured ??.?? μs on the small sample set and ??.?? μs on the large one.
/// </summary>
[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).
/// </summary>
[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;
}
/// <summary>
/// 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).
/// </summary>
[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;
EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs
+77 -8
View File
@@ -89,14 +89,32 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
/// large set.
/// </summary>
// [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;
}
/// <summary>
/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
/// large set.
/// </summary>
// [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×).
/// </summary>
[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;
}
/// <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
/// on the large one (0.08×).
/// </summary>
[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).
/// </summary>
// [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;
}
/// <summary>
/// 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).
/// </summary>
// [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;
EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs
+3 -3
View File
@@ -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)
{
PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+12 -15
View File
@@ -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<DivisorScanSession> _sessionPool = new();
private readonly ConcurrentBag<DivisorScanSession> _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<ulong, MersenneDivisorCycles.FactorCacheEntry>(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)
PerfectNumbers.Core/ExponentRemainderStepper.cs
+16 -27
View File
@@ -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;
PerfectNumbers.Core/FactorEntry.cs
+8
View File
@@ -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;
}
}
PerfectNumbers.Core/Gpu/GpuUInt128.cs
+12 -12
View File
@@ -267,7 +267,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
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<GpuUInt128>, IEquatable<GpuUInt128>
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<GpuUInt128>, IEquatable<GpuUInt128>
}
[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<GpuUInt128>, IEquatable<GpuUInt128>
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<GpuUInt128>, IEquatable<GpuUInt128>
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<GpuUInt128>, IEquatable<GpuUInt128>
}
[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<GpuUInt128>, IEquatable<GpuUInt128>
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<GpuUInt128>, IEquatable<GpuUInt128>
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Mul(in GpuUInt128 other)
public void Mul(GpuUInt128 other)
{
ReadOnlyGpuUInt128 readOnlyThis = AsReadOnly();
ReadOnlyGpuUInt128 readOnlyOther = other.AsReadOnly();
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+2 -2
View File
@@ -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
PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
+4 -4
View File
@@ -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<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1811,7 +1811,7 @@ internal static partial class PrimeOrderGpuHeuristics
private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
{
ulong exponent = exponents[index];
remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
}
private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)
PerfectNumbers.Core/HeuristicFailureLog.cs
+98
View File
@@ -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<StringBuilder?> s_threadBuilder = new(() => StringBuilderPool.Rent(), trackAllValues: true);
private static readonly ThreadLocal<int> 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);
}
}
}
}
PerfectNumbers.Core/HeuristicFailureReason.cs
+11
View File
@@ -0,0 +1,11 @@
namespace PerfectNumbers.Core
{
internal enum HeuristicFailureReason
{
PhiPartialFactorizationFailed,
StrictModeRequested,
HeuristicPipelineUnresolved,
HeuristicPipelineAssumedFullOrder,
}
}
PerfectNumbers.Core/MersenneDivisorCycles.cs
+12 -12
View File
@@ -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<ulong, FactorCacheEntry>? 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<ulong, FactorCacheEntry>? factorCache)
public static bool CycleEqualsExponentCpu(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary<ulong, FactorCacheEntry>? 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<ulong, int> factorCounts)
private static ulong ReduceOrderCpu(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> 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;
PerfectNumbers.Core/PartialFactorResult.cs
+76
View File
@@ -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<FactorEntry>.Shared.Rent(1);
local[0] = new FactorEntry(prime, 1);
return Rent(local, 1UL, true, 1);
}
FactorEntry[] extended = ArrayPool<FactorEntry>.Shared.Rent(Count + 1);
Array.Copy(source, 0, extended, 0, Count);
extended[Count] = new FactorEntry(prime, 1);
Span<FactorEntry> 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<FactorEntry>.Shared.Return(factors, clearArray: false);
}
_next = s_poolHead;
s_poolHead = this;
}
}
}
}
PerfectNumbers.Core/PerfectNumberConstants.cs
+1
View File
@@ -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
PerfectNumbers.Core/PrimeOrderCalculator.cs
+71 -291
View File
@@ -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<ulong, int>? 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<FactorEntry> 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<StringBuilder> s_builderPool = new();
private static readonly ThreadLocal<StringBuilder?> s_threadBuilder = new(() => AcquireBuilder(), trackAllValues: true);
private static readonly ThreadLocal<int> 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<ulong> AcquirePendingCompositeList(int capacityHint)
@@ -952,16 +825,15 @@ internal static partial class PrimeOrderCalculator
private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
{
Dictionary<ulong, int>? dictionary = s_partialFactorCounts;
int adjustedCapacity = Math.Max(capacityHint, 8);
if (dictionary is null)
{
dictionary = new Dictionary<ulong, int>(adjustedCapacity);
dictionary = new Dictionary<ulong, int>(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<ulong, int> 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<FactorEntry>.Shared.Rent(1);
local[0] = new FactorEntry(prime, 1);
return Rent(local, 1UL, true, 1);
}
FactorEntry[] extended = ArrayPool<FactorEntry>.Shared.Rent(Count + 1);
Array.Copy(source, 0, extended, 0, Count);
extended[Count] = new FactorEntry(prime, 1);
Span<FactorEntry> 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<FactorEntry>.Shared.Return(factors, clearArray: false);
}
Cofactor = 0UL;
FullyFactored = false;
Count = 0;
_next = s_poolHead;
s_poolHead = this;
}
}
}
PerfectNumbers.Core/UInt128Extensions.cs
+6 -8
View File
@@ -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)]
PerfectNumbers.Core/ULongExtensions.cs
+117 -8
View File
@@ -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<ulong>.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<ulong>.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)
{