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

Remove default prime order overloads

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-10-10 21:51:41 +02:00
1 parent 64c08f9b49
commit c9ae56a5ec
28 files changed
+1438 -964
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+14 -15
View File
@@ -238,20 +238,19 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
}
if (divisorCycle != 0 && divisorCycle == prime)
{
processedAll = true;
return true;
}
else if (divisorCycle == 0)
{
Console.WriteLine($"Divisor cycle wasn't calculated for ${prime}");
}
// if (divisorCycle != 0UL && CheckDivisor(prime, divisorCycle, divisorData) != 0)
// {
// processedAll = true;
// return true;
// }
if (divisorCycle == prime)
{
// A cycle equal to the tested exponent (which is prime in this path)
// guarantees that the candidate divides the corresponding Mersenne
// number because the order of 2 modulo the divisor is exactly p.
processedAll = true;
return true;
}
if (divisorCycle == 0UL)
{
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
}
}
divisor += step;
@@ -269,7 +268,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;
}
PerfectNumbers.Core/DeterministicRandom.cs
+49
View File
@@ -0,0 +1,49 @@
using System;
using System.Threading;
namespace PerfectNumbers.Core;
internal static class DeterministicRandom
{
[ThreadStatic]
private static ulong s_threadState;
private static long s_seedCounter;
public static ulong NextUInt64()
{
ulong state = s_threadState;
if (state == 0UL)
{
state = InitializeThreadState();
}
state ^= state >> 12;
state ^= state << 25;
state ^= state >> 27;
s_threadState = state;
return state * 2685821657736338717UL;
}
public static UInt128 NextUInt128()
{
UInt128 high = (UInt128)NextUInt64();
UInt128 low = (UInt128)NextUInt64();
return (high << 64) | low;
}
private static ulong InitializeThreadState()
{
ulong seed = (ulong)Interlocked.Increment(ref s_seedCounter);
seed += 0x9E3779B97F4A7C15UL;
seed = (seed ^ (seed >> 30)) * 0xBF58476D1CE4E5B9UL;
seed = (seed ^ (seed >> 27)) * 0x94D049BB133111EBUL;
seed ^= seed >> 31;
if (seed == 0UL)
{
seed = 0x9E3779B97F4A7C15UL;
}
return seed;
}
}
PerfectNumbers.Core/ExponentRemainderStepper.cs
+5 -5
View File
@@ -47,7 +47,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 +59,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 +77,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 +88,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 +123,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/GpuContextPool.cs
+1 -1
View File
@@ -126,7 +126,7 @@ public static class GpuContextPool
}
}
public struct GpuContextLease : IDisposable
public struct GpuContextLease
{
private readonly PooledContext _ctx;
private bool disposed;
PerfectNumbers.Core/Gpu/GpuKernelPool.cs
+1 -1
View File
@@ -677,7 +677,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
// GC.SuppressFinalize(this);
}
public readonly struct ExecutionScope : IDisposable
public readonly struct ExecutionScope
{
private readonly object? _lock;
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
+11 -12
View File
@@ -415,9 +415,9 @@ internal static partial class PrimeOrderGpuHeuristics
ulong? previousOrder,
PrimeOrderCalculator.PrimeOrderSearchConfig config,
in MontgomeryDivisorData divisorData,
out PrimeOrderCalculator.PrimeOrderResult result)
out ulong order)
{
result = default;
order = 0UL;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
@@ -497,18 +497,17 @@ internal static partial class PrimeOrderGpuHeuristics
throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity.");
}
ulong order = 0UL;
resultBuffer.View.CopyToCPU(ref order, 1);
PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found
? PrimeOrderCalculator.PrimeOrderStatus.Found
: PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved;
if (kernelStatus == PrimeOrderKernelStatus.FactoringFailure)
{
order = 0UL;
}
result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order);
DisposeResources();
lease.Dispose();
return true;
return order != 0UL;
void DisposeResources()
{
@@ -812,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;
}
@@ -1596,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)
@@ -1796,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)
@@ -1812,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/MersenneDivisorCycles.cs
+12 -12
View File
@@ -102,7 +102,7 @@ public class MersenneDivisorCycles
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)
@@ -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)
@@ -234,14 +234,14 @@ public class MersenneDivisorCycles
return (UInt128)cycle;
}
PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
UInt128 wideOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (wideResult.Order != UInt128.Zero)
if (wideOrder != UInt128.Zero)
{
return wideResult.Order;
return wideOrder;
}
// Otherwise, find order of 2 mod divisor
@@ -285,15 +285,15 @@ public class MersenneDivisorCycles
return true;
}
PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (orderResult.Order != 0UL)
if (computedOrder != 0UL)
{
cycleLength = orderResult.Order;
cycleLength = computedOrder;
return true;
}
@@ -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;
@@ -940,15 +940,15 @@ public class MersenneDivisorCycles
if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
{
PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (orderResult.Order != 0UL)
if (computedOrder != 0UL)
{
cycleLength = orderResult.Order;
cycleLength = computedOrder;
return true;
}
}
PerfectNumbers.Core/MersennePrimeFactorTester.cs
+4 -11
View File
@@ -371,10 +371,6 @@ public static class MersennePrimeFactorTester
return 2UL;
}
var rng = new Random(1234567); // TODO: Replace this RNG with the deterministic span-based
// sequence measured faster in PollardRhoBenchmarks so we avoid
// per-call Random allocations.
ulong c, d, x, y, diff;
while (true)
{
@@ -383,8 +379,8 @@ public static class MersennePrimeFactorTester
return 0UL;
}
c = (ulong)rng.Next(2, int.MaxValue);
x = (ulong)rng.Next(2, int.MaxValue);
c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
y = x;
d = 1UL;
@@ -523,9 +519,6 @@ public static class MersennePrimeFactorTester
return 2;
}
var rng = new Random(1234567); // TODO: Replace with the deterministic UInt128 sequence from
// PollardRho128Benchmarks so we eliminate Random allocations.
UInt128 c, d, x, y, diff;
while (true)
{
@@ -534,8 +527,8 @@ public static class MersennePrimeFactorTester
return 0;
}
c = (UInt128)rng.NextInt64(2, long.MaxValue);
x = (UInt128)rng.NextInt64(2, long.MaxValue);
c = (DeterministicRandom.NextUInt128() % (n - 1)) + 1;
x = (DeterministicRandom.NextUInt128() % (n - 2)) + 2;
y = x;
d = 1;
PerfectNumbers.Core/PartialFactorResult.cs
+81
View File
@@ -0,0 +1,81 @@
using System;
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 PartialFactorResult();
}
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 PartialFactorResult
{
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.Wide.cs
+28 -37
View File
@@ -9,30 +9,15 @@ namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
internal readonly struct PrimeOrderResultWide
{
public PrimeOrderResultWide(PrimeOrderStatus status, UInt128 order)
{
Status = status;
Order = order;
}
public PrimeOrderStatus Status { get; }
public UInt128 Order { get; }
}
public static PrimeOrderResultWide Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
=> Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
public static PrimeOrderResultWide Calculate(
public static UInt128 Calculate(
UInt128 prime,
UInt128? previousOrder,
PrimeOrderSearchConfig config,
PrimeOrderHeuristicDevice device)
{
using var scope = UsePow2Mode(device);
var scope = UsePow2Mode(device);
MontgomeryDivisorData divisorData;
UInt128 result;
if (prime <= ulong.MaxValue)
{
ulong? previous = null;
@@ -51,41 +36,49 @@ internal static partial class PrimeOrderCalculator
ulong prime64 = (ulong)prime;
divisorData = MontgomeryDivisorData.FromModulus(prime64);
PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
return new PrimeOrderResultWide(result.Status, result.Order);
ulong order64 = Calculate(prime64, previous, divisorData, config, device);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
}
else
{
divisorData = default;
result = CalculateWideInternal(prime, previousOrder, divisorData, config);
}
return CalculateWideInternal(prime, previousOrder, divisorData, config);
scope.Dispose();
return result;
}
private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
{
if (prime <= UInt128.One)
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, UInt128.One);
return UInt128.One;
}
if (prime == (UInt128)3UL)
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, (UInt128)2UL);
return (UInt128)2UL;
}
UInt128 phi = prime - UInt128.One;
using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
UInt128 result;
if (phiFactors.Factors is null)
{
return FinishStrictlyWide(prime, divisorData, config.Mode);
result = FinishStrictlyWide(prime, divisorData);
}
else
{
result = RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
}
return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
debugScope.Dispose();
return result;
}
private static PrimeOrderResultWide RunHeuristicPipelineWide(
private static UInt128 RunHeuristicPipelineWide(
UInt128 prime,
UInt128? previousOrder,
in MontgomeryDivisorData divisorData,
@@ -95,7 +88,7 @@ internal static partial class PrimeOrderCalculator
{
if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
return phi;
}
UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
@@ -103,27 +96,25 @@ internal static partial class PrimeOrderCalculator
if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder);
return candidateOrder;
}
if (config.Mode == PrimeOrderMode.Strict)
{
return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict);
return FinishStrictlyWide(prime, divisorData);
}
if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, order);
return order;
}
return FinishStrictlyWide(prime, divisorData, config.Mode);
return candidateOrder;
}
private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode)
private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData)
{
UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData);
PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
return new PrimeOrderResultWide(status, strictOrder);
return CalculateByFactorizationWide(prime, divisorData);
}
private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
PerfectNumbers.Core/PrimeOrderCalculator.cs
+545 -468
View File
File diff suppressed because it is too large. Load diff
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 Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
{
return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
{
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)
{