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

Revert "Reusable buffers introduced."

Browse Source 
This reverts commit 39a66fdef1.
This commit is contained in:
mlemanczyk committed 2025-12-13 04:17:30 +01:00
1 parent 39a66fdef1
commit 8fbe248a49
8 files changed
+1056 -829
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
+1 -1
View File
@@ -7,7 +7,7 @@ namespace PerfectNumbers.Core;
/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path.
/// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order.
/// </summary>
public struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL)
internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL)
{
private readonly MontgomeryDivisorData _divisor = divisor;
private readonly ulong _modulus = divisor.Modulus;
PerfectNumbers.Core/PartialFactorPendingEntry.cs
-47
View File
@@ -1,47 +0,0 @@
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
internal sealed class PartialFactorPendingEntry(ulong value, bool knownComposite)
{
[ThreadStatic]
private static PartialFactorPendingEntry? _poolHead;
private PartialFactorPendingEntry? _next;
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static PartialFactorPendingEntry Rent(ulong value, bool knownComposite)
{
var poolHead = _poolHead;
if (poolHead == null)
{
return new(value, knownComposite);
}
_poolHead = poolHead._next;
poolHead.Value = value;
poolHead.KnownComposite = knownComposite;
poolHead.HasKnownPrimality = knownComposite;
poolHead.IsPrime = false;
return poolHead;
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static void Return(PartialFactorPendingEntry entry)
{
entry._next = _poolHead;
_poolHead = entry;
}
public ulong Value = value;
public bool KnownComposite = knownComposite;
public bool HasKnownPrimality = knownComposite;
public bool IsPrime = false;
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public void WithPrimality(bool isPrime)
{
KnownComposite = KnownComposite || !isPrime;
HasKnownPrimality = true;
IsPrime = isPrime;
}
}
PerfectNumbers.Core/PartialFactorResult.cs
+109 -64
View File
@@ -1,109 +1,154 @@
using System.Collections.Generic;
using System.Numerics;
using System.Runtime.CompilerServices;
using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core
{
public sealed class PartialFactorResult
{
private const int ExponentHardLimit = 256;
[ThreadStatic]
private static PartialFactorResult? s_poolHead;
public static readonly PartialFactorResult Empty = new(MontgomeryDivisorData.Empty)
public ulong[]? Factors;
public int[]? Exponents;
private PartialFactorResult? _next;
public ulong Cofactor;
public bool FullyFactored;
public int Count;
public bool CofactorIsPrime;
private PartialFactorResult()
{
}
public static PartialFactorResult Rent(in ulong[]? factors, in int[]? exponents, ulong cofactor, bool fullyFactored, int count, bool cofactorIsPrime)
{
PartialFactorResult? instance = s_poolHead;
if (instance is null)
{
instance = new PartialFactorResult();
}
else
{
s_poolHead = instance._next;
}
instance._next = null;
instance.Factors = factors;
instance.Exponents = exponents;
instance.Cofactor = cofactor;
instance.FullyFactored = fullyFactored;
instance.Count = count;
instance.CofactorIsPrime = cofactorIsPrime;
return instance;
}
public static PartialFactorResult Rent(ulong cofactor, bool fullyFactored, int count, bool cofactorIsPrime)
{
PartialFactorResult? instance = s_poolHead;
if (instance is null)
{
instance = new PartialFactorResult();
}
else
{
s_poolHead = instance._next;
}
instance._next = null;
instance.Factors = null;
instance.Exponents = null;
instance.Cofactor = cofactor;
instance.FullyFactored = fullyFactored;
instance.Count = count;
instance.CofactorIsPrime = cofactorIsPrime;
return instance;
}
public static readonly PartialFactorResult Empty = new()
{
Cofactor = 1UL,
FullyFactored = true
};
internal readonly List<PartialFactorPendingEntry> Pending = new(PrimeOrderConstants.GpuSmallPrimeFactorSlots);
internal readonly FixedCapacityStack<ulong> CompositeStack = new(PrimeOrderConstants.GpuSmallPrimeFactorSlots);
internal readonly ulong[] SpecialMaxBuffer = new ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots];
public readonly ulong[] Factors = new ulong[ExponentHardLimit];
public readonly int[] Exponents = new int[ExponentHardLimit];
public readonly ulong[] FactorCandidates = new ulong[ExponentHardLimit];
public readonly List<ulong> FactorCandidatesList = new(ExponentHardLimit);
public readonly ulong[] ValidateOrderFactorCandidates = new ulong[ExponentHardLimit];
public readonly bool[] FactorCandidateEvaluations = new bool[ExponentHardLimit];
internal readonly Dictionary<ulong, int> FactorCounts = new(capacity: 8);
public ExponentRemainderStepperCpu ExponentRemainderStepper;
private PartialFactorResult? _next;
public ulong Cofactor;
public bool CofactorIsPrime;
public int Count;
public bool FullyFactored;
public bool HasFactors;
private PartialFactorResult(in MontgomeryDivisorData divisorData)
{
ExponentRemainderStepper = new(divisorData);
}
public static PartialFactorResult Rent(in MontgomeryDivisorData divisorData)
{
if (s_poolHead is { } instance)
{
s_poolHead = instance._next;
instance._next = null;
instance.Cofactor = 0UL;
instance.FullyFactored = false;
instance.Count = 0;
instance.CofactorIsPrime = false;
instance.HasFactors = false;
if (!instance.ExponentRemainderStepper.MatchesDivisor(divisorData, 0UL))
{
instance.ExponentRemainderStepper = new ExponentRemainderStepperCpu(divisorData, 0UL);
}
return instance;
}
return new PartialFactorResult(divisorData);
}
public void WithAdditionalPrime(ulong prime)
{
ulong[] sourceFactors = Factors;
ulong[]? sourceFactors = Factors;
int sourceCount = Count;
if (sourceCount == 0)
if (sourceFactors is null || sourceCount == 0)
{
InitializeWithPrime(prime);
return;
}
int[] sourceExponents = Exponents;
int newSize = (int)BitOperations.RoundUpToPowerOf2((uint)sourceCount + 1U);
ulong[] extendedFactors = sourceFactors.Length >= newSize
? sourceFactors
: ResizeFactors(sourceFactors, sourceCount, newSize);
int newSize = sourceCount + 1;
sourceFactors[sourceCount] = prime;
sourceExponents[sourceCount] = 1;
Array.Sort(sourceFactors, sourceExponents, 0, newSize);
int[] sourceExponents = Exponents!;
int[] extendedExponents = sourceExponents.Length >= newSize
? sourceExponents
: ResizeExponents(sourceCount, newSize, sourceExponents);
newSize = sourceCount + 1;
extendedFactors[sourceCount] = prime;
extendedExponents[sourceCount] = 1;
Array.Sort(extendedFactors, extendedExponents, 0, newSize);
Factors = extendedFactors;
Exponents = extendedExponents;
Cofactor = 1UL;
FullyFactored = true;
Count = newSize;
CofactorIsPrime = false;
HasFactors = true;
}
private static int[] ResizeExponents(int sourceCount, int newSize, in int[] sourceExponents)
{
int[] extendedExponents = FixedCapacityPools.ExclusiveIntArray.Rent(newSize);
Array.Copy(sourceExponents, 0, extendedExponents, 0, sourceCount);
FixedCapacityPools.ExclusiveIntArray.Return(sourceExponents);
return extendedExponents;
}
private static ulong[] ResizeFactors(in ulong[] sourceFactors, int sourceCount, int newSize)
{
ulong[] extendedFactors = FixedCapacityPools.ExclusiveUlongArray.Rent(newSize);
Array.Copy(sourceFactors, 0, extendedFactors, 0, sourceCount);
FixedCapacityPools.ExclusiveUlongArray.Return(sourceFactors);
return extendedFactors;
}
private void InitializeWithPrime(ulong prime)
{
Factors[0] = prime;
Exponents[0] = 1;
ulong[] factors = FixedCapacityPools.ExclusiveUlongArray.Rent(1);
int[] exponents = FixedCapacityPools.ExclusiveIntArray.Rent(1);
factors[0] = prime;
exponents[0] = 1;
Factors = factors;
Exponents = exponents;
Cofactor = 1UL;
FullyFactored = true;
Count = 1;
CofactorIsPrime = false;
HasFactors = true;
}
public void Dispose()
{
if (this != Empty)
{
ulong[]? factors = Factors;
if (factors is not null)
{
int[] exponents = Exponents!;
Factors = null;
Exponents = null;
FixedCapacityPools.ExclusiveUlongArray.Return(factors);
FixedCapacityPools.ExclusiveIntArray.Return(exponents);
}
_next = s_poolHead;
s_poolHead = this;
}
PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
+445 -283
View File
File diff suppressed because it is too large. Load diff
PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs
+177 -147
View File
@@ -76,12 +76,12 @@ internal static partial class PrimeOrderCalculator
return gpuOrder;
}
PartialFactorResult phiFactors = PartialFactorGpu(gpu, divisorData, phi, config);
PartialFactorResult phiFactors = PartialFactorGpu(gpu, phi, config);
ulong result;
if (phiFactors.Count == 0)
if (phiFactors.Factors is null)
{
result = CalculateByFactorizationCpu(prime, divisorData, phiFactors);
result = CalculateByFactorizationCpu(prime, divisorData);
phiFactors.Dispose();
return result;
@@ -147,28 +147,21 @@ internal static partial class PrimeOrderCalculator
UInt128 phi = prime - UInt128.One;
PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
try
UInt128 result;
if (phiFactors.Factors is null)
{
UInt128 result;
if (phiFactors.Count == 0)
{
result = FinishStrictlyWideGpu(gpu, prime, divisorData);
}
else
{
result = RunHeuristicPipelineWideGpu(gpu, prime, previousOrder, divisorData, config, phi, phiFactors);
}
result = FinishStrictlyWideGpu(gpu, prime, divisorData);
}
else
{
result = RunHeuristicPipelineWideGpu(gpu, prime, previousOrder, divisorData, config, phi, phiFactors);
}
return result;
}
finally
{
phiFactors.Dispose();
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static void CollectFactorsGpu(PrimeOrderCalculatorAccelerator gpu, ref Span<ulong> primeSlots, ref Span<int> exponentSlots, ref int factorCount, List<PartialFactorPendingEntry> pending, Stack<ulong> compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached)
private static void CollectFactorsGpu(PrimeOrderCalculatorAccelerator gpu, ref Span<ulong> primeSlots, ref Span<int> exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List<PendingEntry> pending, Stack<ulong> compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached)
{
while (compositeStack.Count > 0)
{
@@ -177,7 +170,7 @@ internal static partial class PrimeOrderCalculator
if (isPrime)
{
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite);
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1);
continue;
}
@@ -185,14 +178,14 @@ internal static partial class PrimeOrderCalculator
if (currentTimestamp > deadlineTimestamp)
{
pollardRhoDeadlineReached = true;
pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true));
pending.Add(new PendingEntry(composite, knownComposite: true));
continue;
}
// if (!TryPollardRhoGpu(gpu, composite, out ulong factor))
if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor))
{
pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true));
pending.Add(new PendingEntry(composite, knownComposite: true));
continue;
}
@@ -207,6 +200,7 @@ internal static partial class PrimeOrderCalculator
PrimeOrderCalculatorAccelerator gpu,
ReadOnlySpan<ulong> factors,
ulong phi,
ulong prime,
in MontgomeryDivisorData divisorData)
{
int factorCount = factors.Length;
@@ -241,6 +235,7 @@ internal static partial class PrimeOrderCalculator
private static UInt128 ExponentLoweringWideGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors)
{
ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool;
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
int length = factors.Count;
// TODO: This should never trigger from production code - check
@@ -248,7 +243,9 @@ internal static partial class PrimeOrderCalculator
// {
// return order;
// }
Span<FactorEntry128> buffer = factors.TempFactors.AsSpan(0, length + 1);
FactorEntry128[] tempArray = pool.Rent(length + 1);
Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
factorSpan.CopyTo(buffer);
if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor))
@@ -266,7 +263,7 @@ internal static partial class PrimeOrderCalculator
int exponent = buffer[i].Exponent;
for (int iteration = 0; iteration < exponent; iteration++)
{
if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero)
if ((order % primeFactor) == UInt128.Zero)
{
UInt128 reduced = order / primeFactor;
if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One)
@@ -280,6 +277,7 @@ internal static partial class PrimeOrderCalculator
}
}
pool.Return(tempArray, clearArray: false);
return order;
}
@@ -304,7 +302,7 @@ internal static partial class PrimeOrderCalculator
int exponent = entries[i].Value;
for (int iteration = 0; iteration < exponent; iteration++)
{
if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero)
if ((order % primeFactor) != UInt128.Zero)
{
UInt128 candidate = order / primeFactor;
if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) == UInt128.One)
@@ -339,26 +337,28 @@ internal static partial class PrimeOrderCalculator
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisorData, ulong value, in PrimeOrderCalculatorConfig config)
private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config)
{
if (value <= 1UL)
{
{
return PartialFactorResult.Empty;
}
// stackalloc is faster than pooling
PartialFactorResult result = PartialFactorResult.Rent(divisorData);
Span<ulong> primeSlots = new(result.Factors);
Span<int> exponentSlots = new (result.Exponents);
Span<ulong> primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots];
Span<int> exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots];
ulong[] primeSlotArray = [];
int[] exponentSlotArray = [];
// We don't need to worry about leftovers, because we always use indexes within the calculated counts
// primeSlots.Clear();
// exponentSlots.Clear();
uint limit = config.SmallFactorLimit;
bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out ulong remaining);
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
ulong remaining = value;
List<PartialFactorPendingEntry> pending = result.Pending;
pending.Clear();
bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out remaining);
List<PendingEntry> pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2);
if (!gpuFactored)
{
@@ -371,28 +371,31 @@ internal static partial class PrimeOrderCalculator
out remaining);
}
bool limitReached;
if (remaining > 1UL)
{
if (config.PollardRhoMilliseconds > 0)
{
long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
FixedCapacityStack<ulong> compositeStack = result.CompositeStack;
compositeStack.Clear();
FixedCapacityStack<ulong> compositeStack = ThreadStaticPools.RentUlongStack(4);
compositeStack.Push(remaining);
CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out bool limitReached);
CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached);
if (limitReached)
{
while (compositeStack.Count > 0)
{
pending.Add(PartialFactorPendingEntry.Rent(compositeStack.Pop(), knownComposite: false));
pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false));
}
}
compositeStack.Clear();
ThreadStaticPools.ReturnUlongStack(compositeStack);
}
else
{
pending.Add(PartialFactorPendingEntry.Rent(remaining, knownComposite: false));
pending.Add(new PendingEntry(remaining, knownComposite: false));
}
}
@@ -402,14 +405,13 @@ internal static partial class PrimeOrderCalculator
int index = 0;
for (; index < pendingCount; index++)
{
PartialFactorPendingEntry entry = pending[index];
PendingEntry entry = pending[index];
ulong composite = entry.Value;
if (entry.KnownComposite)
{
cofactor = checked(cofactor * composite);
cofactorContainsComposite = true;
PartialFactorPendingEntry.Return(entry);
continue;
}
@@ -423,7 +425,7 @@ internal static partial class PrimeOrderCalculator
if (entry.IsPrime)
{
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite);
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1);
}
// composite is never smaller on the execution path
// else if (composite > 1UL)
@@ -432,8 +434,6 @@ internal static partial class PrimeOrderCalculator
cofactor = checked(cofactor * composite);
cofactorContainsComposite = true;
}
PartialFactorPendingEntry.Return(entry);
}
bool cofactorIsPrime;
@@ -450,6 +450,19 @@ internal static partial class PrimeOrderCalculator
cofactorIsPrime = PrimeTesterByLastDigit.IsPrimeCpu(cofactor);
}
PartialFactorResult result;
if (factorCount == 0)
{
if (cofactor == value)
{
result = PartialFactorResult.Rent(cofactor, false, 0, cofactorIsPrime);
goto ReturnResult;
}
result = PartialFactorResult.Rent(cofactor, cofactor == 1UL, 0, cofactorIsPrime);
goto ReturnResult;
}
// This will never happen in production code. We'll always get at least 1 factor
// if (factorCount == 0)
// {
@@ -463,13 +476,41 @@ internal static partial class PrimeOrderCalculator
// goto ReturnResult;
// }
Array.Sort(result.Factors, result.Exponents, 0, factorCount);
ArrayPool<ulong> ulongPool = ThreadStaticPools.UlongPool;
ArrayPool<int> intPool = ThreadStaticPools.IntPool;
ulong[] factorsArray = ulongPool.Rent(factorCount);
int[] exponentsArray = intPool.Rent(factorCount);
for (int i = 0; i < factorCount; i++)
{
factorsArray[i] = primeSlots[i];
exponentsArray[i] = exponentSlots[i];
// This will never happen on the execution path from production code
// if (primeValue == 0UL || exponentValue == 0)
// {
// throw new Exception("Prime value or exponent equals zero");
// continue;
// }
}
Array.Sort(factorsArray, exponentsArray, 0, factorCount);
// Check if it was factored completely
result.Cofactor = cofactor;
result.FullyFactored = cofactor == 1UL;
result.Count = factorCount;
result.CofactorIsPrime = cofactorIsPrime;
limitReached = cofactor == 1UL;
result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime);
ReturnResult:
if (primeSlotArray.Length > 0)
{
ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false);
ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false);
}
pending.Clear();
ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending);
return result;
}
@@ -539,9 +580,8 @@ internal static partial class PrimeOrderCalculator
private static bool TryConfirmCandidateWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData? divisorData)
{
PartialFactorResult128 factorization = PartialFactorWide(candidate, config);
if (factorization.Count == 0)
if (factorization.Factors is null)
{
factorization.Dispose();
return false;
}
@@ -549,54 +589,45 @@ internal static partial class PrimeOrderCalculator
{
if (factorization.Cofactor <= UInt128.One)
{
factorization.Dispose();
return false;
}
if (!IsPrimeWide(factorization.Cofactor))
{
factorization.Dispose();
return false;
}
factorization.WithAdditionalPrime(factorization.Cofactor);
factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
}
ReadOnlySpan<FactorEntry128> span = factorization.Factors;
int length = factorization.Count;
try
for (int i = 0; i < length; i++)
{
for (int i = 0; i < length; i++)
UInt128 primeFactor = span[i].Value;
UInt128 reduced = candidate;
for (int iteration = 0; iteration < span[i].Exponent; iteration++)
{
UInt128 primeFactor = span[i].Value;
UInt128 reduced = candidate;
for (int iteration = 0; iteration < span[i].Exponent; iteration++)
if ((reduced % primeFactor) != UInt128.Zero)
{
if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero)
{
break;
}
break;
}
reduced /= primeFactor;
if (powUsed >= powBudget && powBudget > 0)
{
return false;
}
reduced /= primeFactor;
if (powUsed >= powBudget && powBudget > 0)
{
return false;
}
powUsed++;
if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One)
{
return false;
}
powUsed++;
if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One)
{
return false;
}
}
}
return true;
}
finally
{
factorization.Dispose();
}
return true;
}
private static bool TryConfirmOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config)
@@ -612,9 +643,8 @@ internal static partial class PrimeOrderCalculator
}
PartialFactorResult128 factorization = PartialFactorWide(order, config);
if (factorization.Count == 0)
if (factorization.Factors is null)
{
factorization.Dispose();
return false;
}
@@ -622,17 +652,15 @@ internal static partial class PrimeOrderCalculator
{
if (factorization.Cofactor <= UInt128.One)
{
factorization.Dispose();
return false;
}
if (!IsPrimeWide(factorization.Cofactor))
{
factorization.Dispose();
return false;
}
factorization.WithAdditionalPrime(factorization.Cofactor);
factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
}
ReadOnlySpan<FactorEntry128> span = factorization.Factors;
@@ -651,13 +679,11 @@ internal static partial class PrimeOrderCalculator
reduced /= primeFactor;
if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One)
{
factorization.Dispose();
return false;
}
}
}
factorization.Dispose();
return true;
}
@@ -669,6 +695,7 @@ internal static partial class PrimeOrderCalculator
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config,
PartialFactorResult phiFactors,
PartialFactorResult? cachedOrderFactors,
out ulong result)
{
@@ -680,10 +707,10 @@ internal static partial class PrimeOrderCalculator
}
// Reuse the partial factorization from TryConfirmOrderCpu when available.
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, divisorData, order, config);
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, config);
try
{
if (orderFactors.Count == 0)
if (orderFactors.Factors is null)
{
return false;
}
@@ -699,13 +726,16 @@ internal static partial class PrimeOrderCalculator
}
int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity;
List<ulong> candidates = orderFactors.FactorCandidatesList;
List<ulong> candidates = ThreadStaticPools.RentUlongList(capacity);
candidates.Clear();
ulong[] factorArray = orderFactors.Factors!;
int[] exponentsArray = orderFactors.Exponents!;
BuildCandidates(order, orderFactors, orderFactors.Count, candidates, capacity);
BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity);
int candidateCount = candidates.Count;
if (candidateCount == 0)
{
ThreadStaticPools.ReturnUlongList(candidates);
return false;
}
@@ -715,12 +745,11 @@ internal static partial class PrimeOrderCalculator
int powUsed = 0;
bool allowGpuBatch = true;
Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
ExponentRemainderStepperCpu powStepper = orderFactors.ExponentRemainderStepper;
powStepper.Reset();
ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool powStepperInitialized = false;
int index = 0;
Span<ulong> stackGpuRemainders = new(orderFactors.SpecialMaxBuffer, 0, PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount);
Span<ulong> stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount];
ArrayPool<ulong> pool = ThreadStaticPools.UlongPool;
while (index < candidateCount && powUsed < powBudget)
{
@@ -790,6 +819,9 @@ internal static partial class PrimeOrderCalculator
pool.Return(gpuPool, clearArray: false);
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
@@ -827,6 +859,9 @@ internal static partial class PrimeOrderCalculator
pool.Return(gpuPool, clearArray: false);
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
@@ -839,6 +874,9 @@ internal static partial class PrimeOrderCalculator
index += batchSize;
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
// DebugLog("No candidate confirmed");
return false;
}
@@ -864,75 +902,66 @@ internal static partial class PrimeOrderCalculator
}
PartialFactorResult128 orderFactors = PartialFactorWide(order, config);
if (orderFactors.Count == 0)
if (orderFactors.Factors is null)
{
orderFactors.Dispose();
return false;
}
try
if (!orderFactors.FullyFactored)
{
if (!orderFactors.FullyFactored)
{
if (orderFactors.Cofactor <= UInt128.One)
{
return false;
}
if (!IsPrimeWide(orderFactors.Cofactor))
{
return false;
}
orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
}
int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity;
List<UInt128> candidates = orderFactors.FactorCandidatesList;
candidates.Clear();
FactorEntry128[] factorArray = orderFactors.Factors!;
BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity);
if (candidates.Count == 0)
if (orderFactors.Cofactor <= UInt128.One)
{
return false;
}
SortCandidatesWide(prime, previousOrder, candidates);
int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
int powUsed = 0;
int candidateCount = candidates.Count;
for (int i = 0; i < candidateCount; i++)
if (!IsPrimeWide(orderFactors.Cofactor))
{
if (powUsed >= powBudget && powBudget > 0)
{
break;
}
UInt128 candidate = candidates[i];
powUsed++;
if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) != UInt128.One)
{
continue;
}
if (!TryConfirmCandidateWideGpu(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData))
{
continue;
}
result = candidate;
return true;
return false;
}
orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
}
int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity;
List<UInt128> candidates = new(capacity);
FactorEntry128[] factorArray = orderFactors.Factors!;
BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity);
if (candidates.Count == 0)
{
return false;
}
finally
SortCandidatesWide(prime, previousOrder, candidates);
int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
int powUsed = 0;
int candidateCount = candidates.Count;
for (int i = 0; i < candidateCount; i++)
{
orderFactors.Dispose();
if (powUsed >= powBudget && powBudget > 0)
{
break;
}
UInt128 candidate = candidates[i];
powUsed++;
if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) != UInt128.One)
{
continue;
}
if (!TryConfirmCandidateWideGpu(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData))
{
continue;
}
result = candidate;
return true;
}
return false;
}
private static bool TryPollardRhoGpu(PrimeOrderCalculatorAccelerator gpu, ulong n, out ulong factor)
@@ -1048,3 +1077,4 @@ internal static partial class PrimeOrderCalculator
return true;
}
}
PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs
+179 -83
View File
@@ -1,6 +1,8 @@
using System.Buffers;
using System.Diagnostics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using ILGPU;
using PerfectNumbers.Core.Gpu;
using PerfectNumbers.Core.Gpu.Accelerators;
@@ -23,12 +25,12 @@ internal static partial class PrimeOrderCalculator
ulong phi = prime - 1UL;
PartialFactorResult phiFactors = PartialFactorHybrid(gpu, divisorData, phi, config);
PartialFactorResult phiFactors = PartialFactorHybrid(gpu, phi, config);
ulong result;
if (phiFactors.Count == 0)
if (phiFactors.Factors is null)
{
result = CalculateByFactorizationCpu(prime, divisorData, phiFactors);
result = CalculateByFactorizationCpu(prime, divisorData);
phiFactors.Dispose();
return result;
@@ -73,14 +75,15 @@ internal static partial class PrimeOrderCalculator
}
else
{
result = CalculateWideInternalCpu(prime, previousOrder, config);
divisorData = MontgomeryDivisorData.Empty;
result = CalculateWideInternalCpu(prime, previousOrder, divisorData, config);
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span<ulong> primeSlots, Span<int> exponentSlots, ref int factorCount, List<PartialFactorPendingEntry> pending, FixedCapacityStack<ulong> compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached)
private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span<ulong> primeSlots, Span<int> exponentSlots, ref int factorCount, List<PendingEntry> pending, FixedCapacityStack<ulong> compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached)
{
while (compositeStack.Count > 0)
{
@@ -89,20 +92,20 @@ internal static partial class PrimeOrderCalculator
if (isPrime)
{
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite);
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1);
continue;
}
pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp;
if (pollardRhoDeadlineReached)
{
pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true));
pending.Add(new PendingEntry(composite, knownComposite: true));
continue;
}
if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor))
{
pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true));
pending.Add(new PendingEntry(composite, knownComposite: true));
continue;
}
@@ -131,7 +134,7 @@ internal static partial class PrimeOrderCalculator
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisorData, ulong value, in PrimeOrderCalculatorConfig config)
private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config)
{
if (value <= 1UL)
{
@@ -139,17 +142,18 @@ internal static partial class PrimeOrderCalculator
}
// stackalloc is faster than pooling
PartialFactorResult result = PartialFactorResult.Rent(divisorData);
Span<ulong> primeSlots = new (result.Factors);
Span<int> exponentSlots = new (result.Exponents);
Span<ulong> primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots];
Span<int> exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots];
ulong[] primeSlotArray = [];
int[] exponentSlotArray = [];
// We don't need to worry about leftovers, because we always use indexes within the calculated counts
// primeSlots.Clear();
// exponentSlots.Clear();
uint limit = config.SmallFactorLimit;
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
ulong remaining = value;
List<PartialFactorPendingEntry> pending = result.Pending;
pending.Clear();
List<PendingEntry> pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2);
PopulateSmallPrimeFactorsCpu(
value,
@@ -157,30 +161,33 @@ internal static partial class PrimeOrderCalculator
primeSlots,
exponentSlots,
out int factorCount,
out ulong remaining);
out remaining);
bool limitReached;
if (remaining > 1UL)
{
if (config.PollardRhoMilliseconds > 0)
{
long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
FixedCapacityStack<ulong> compositeStack = result.CompositeStack;
compositeStack.Clear();
FixedCapacityStack<ulong> compositeStack = ThreadStaticPools.RentUlongStack(4);
compositeStack.Push(remaining);
CollectFactorsHybrid(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out bool limitReached);
CollectFactorsHybrid(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached);
if (limitReached)
{
while (compositeStack.Count > 0)
{
pending.Add(PartialFactorPendingEntry.Rent(compositeStack.Pop(), knownComposite: false));
pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false));
}
}
compositeStack.Clear();
ThreadStaticPools.ReturnUlongStack(compositeStack);
}
else
{
pending.Add(PartialFactorPendingEntry.Rent(remaining, knownComposite: false));
pending.Add(new PendingEntry(remaining, knownComposite: false));
}
}
@@ -190,14 +197,13 @@ internal static partial class PrimeOrderCalculator
int index = 0;
for (; index < pendingCount; index++)
{
PartialFactorPendingEntry entry = pending[index];
PendingEntry entry = pending[index];
ulong composite = entry.Value;
if (entry.KnownComposite)
{
cofactor = checked(cofactor * composite);
cofactorContainsComposite = true;
PartialFactorPendingEntry.Return(entry);
continue;
}
@@ -211,7 +217,7 @@ internal static partial class PrimeOrderCalculator
if (entry.IsPrime)
{
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite);
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1);
}
// composite is never smaller on the execution path
// else if (composite > 1UL)
@@ -220,8 +226,6 @@ internal static partial class PrimeOrderCalculator
cofactor = checked(cofactor * composite);
cofactorContainsComposite = true;
}
PartialFactorPendingEntry.Return(entry);
}
bool cofactorIsPrime;
@@ -238,16 +242,20 @@ internal static partial class PrimeOrderCalculator
cofactorIsPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(cofactor) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor);
}
// This will never happen in production code. We'll always get at least 1 factor
// if (factorCount == 0)
// {
// result.Cofactor = cofactor;
// result.FullyFactored = cofactor != value && cofactor == 1UL;
// result.Count = 0;
// result.CofactorIsPrime = cofactorIsPrime;
// goto ReturnResult;
// }
PartialFactorResult result;
if (factorCount == 0)
{
if (cofactor == value)
{
result = PartialFactorResult.Rent(cofactor, false, 0, cofactorIsPrime);
goto ReturnResult;
}
result = PartialFactorResult.Rent(cofactor, cofactor == 1UL, 0, cofactorIsPrime);
goto ReturnResult;
}
// This will never happen in production code. We'll always get at least 1 factor
// if (factorCount == 0)
// {
// if (cofactor == value)
@@ -260,13 +268,41 @@ internal static partial class PrimeOrderCalculator
// goto ReturnResult;
// }
Array.Sort(result.Factors, result.Exponents, 0, factorCount);
ArrayPool<ulong> ulongPool = ThreadStaticPools.UlongPool;
ArrayPool<int> intPool = ThreadStaticPools.IntPool;
ulong[] factorsArray = ulongPool.Rent(factorCount);
int[] exponentsArray = intPool.Rent(factorCount);
for (int i = 0; i < factorCount; i++)
{
factorsArray[i] = primeSlots[i];
exponentsArray[i] = exponentSlots[i];
// This will never happen on the execution path from production code
// if (primeValue == 0UL || exponentValue == 0)
// {
// throw new Exception("Prime value or exponent equals zero");
// continue;
// }
}
Array.Sort(factorsArray, exponentsArray, 0, factorCount);
// Check if it was factored completely
result.Cofactor = cofactor;
result.FullyFactored = cofactor == 1UL;
result.Count = factorCount;
result.CofactorIsPrime = cofactorIsPrime;
limitReached = cofactor == 1UL;
result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime);
ReturnResult:
if (primeSlotArray.Length > 0)
{
ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false);
ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false);
}
pending.Clear();
ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending);
return result;
}
@@ -300,13 +336,13 @@ internal static partial class PrimeOrderCalculator
ulong phi,
PartialFactorResult phiFactors)
{
if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, phiFactors, divisorData))
if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, prime, phiFactors, divisorData))
{
return phi;
}
ulong candidateOrder = InitializeStartingOrderHybrid(gpu, prime, phi, divisorData);
candidateOrder = ExponentLoweringCpu(candidateOrder, phiFactors, divisorData);
candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData);
if (TryConfirmOrderHybrid(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors))
{
@@ -316,10 +352,10 @@ internal static partial class PrimeOrderCalculator
if (config.StrictMode)
{
orderFactors?.Dispose();
return CalculateByFactorizationCpu(prime, divisorData, phiFactors);
return CalculateByFactorizationCpu(prime, divisorData);
}
if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, orderFactors, out ulong order))
if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order))
{
return order;
}
@@ -371,9 +407,9 @@ internal static partial class PrimeOrderCalculator
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget)
{
PartialFactorResult factorization = PartialFactorHybrid(gpu, divisorData, candidate, config);
PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, config);
if (factorization.Count == 0)
if (factorization.Factors is null)
{
factorization.Dispose();
return false;
@@ -394,8 +430,12 @@ internal static partial class PrimeOrderCalculator
ReadOnlySpan<int> exponentsSpan = factorization.Exponents;
int length = factorization.Count;
Span<ulong> stackBuffer = new(factorization.SpecialMaxBuffer);
ExponentRemainderStepperCpu stepper = factorization.ExponentRemainderStepper;
const int StackExponentCapacity = 32;
const int ExponentHardLimit = 256;
Span<ulong> stackBuffer = stackalloc ulong[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool violates = false;
for (int i = 0; i < length; i++)
@@ -409,13 +449,35 @@ internal static partial class PrimeOrderCalculator
// continue;
// }
if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
if (exponent > ExponentHardLimit)
{
throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}.");
}
if (exponent <= StackExponentCapacity)
{
if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
{
violates = true;
break;
}
continue;
}
heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit);
if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
{
violates = true;
break;
}
}
ThreadStaticPools.ReturnExponentStepperCpu(stepper);
if (heapCandidateArray is not null)
{
ThreadStaticPools.UlongPool.Return(heapCandidateArray);
}
factorization.Dispose();
return !violates;
@@ -442,8 +504,8 @@ internal static partial class PrimeOrderCalculator
return false;
}
PartialFactorResult factorization = PartialFactorHybrid(gpu, divisorData, order, config);
if (factorization.Count == 0)
PartialFactorResult factorization = PartialFactorHybrid(gpu, order, config);
if (factorization.Factors is null)
{
factorization.Dispose();
return false;
@@ -460,7 +522,7 @@ internal static partial class PrimeOrderCalculator
factorization.WithAdditionalPrime(factorization.Cofactor);
}
if (!ValidateOrderAgainstFactors(order, factorization))
if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization))
{
reusableFactorization = factorization;
return false;
@@ -478,6 +540,7 @@ internal static partial class PrimeOrderCalculator
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config,
PartialFactorResult phiFactors,
PartialFactorResult? cachedOrderFactors,
out ulong result)
{
@@ -489,10 +552,10 @@ internal static partial class PrimeOrderCalculator
}
// Reuse the partial factorization from TryConfirmOrderCpu when available.
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, divisorData, order, config);
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, config);
try
{
if (orderFactors.Count == 0)
if (orderFactors.Factors is null)
{
return false;
}
@@ -508,13 +571,16 @@ internal static partial class PrimeOrderCalculator
}
int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity;
List<ulong> candidates = orderFactors.FactorCandidatesList;
List<ulong> candidates = ThreadStaticPools.RentUlongList(capacity);
candidates.Clear();
ulong[] factorArray = orderFactors.Factors!;
int[] exponentsArray = orderFactors.Exponents!;
BuildCandidates(order, orderFactors, orderFactors.Count, candidates, capacity);
BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity);
int candidateCount = candidates.Count;
if (candidateCount == 0)
{
ThreadStaticPools.ReturnUlongList(candidates);
return false;
}
@@ -523,46 +589,74 @@ internal static partial class PrimeOrderCalculator
int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks;
int powUsed = 0;
Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
ExponentRemainderStepperCpu powStepper = orderFactors.ExponentRemainderStepper;
powStepper.Reset();
ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool powStepperInitialized = false;
ulong candidate = candidateSpan[0];
powUsed++;
bool equalsOne = powStepper.InitializeCpuIsUnity(candidate);
if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget))
{
result = candidate;
return true;
}
int index = 1;
int index = 0;
Span<ulong> stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount];
ArrayPool<ulong> pool = ThreadStaticPools.UlongPool;
while (index < candidateCount && powUsed < powBudget)
{
int remaining = candidateCount - index;
int budgetRemaining = powBudget - powUsed;
if (budgetRemaining <= 0)
int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize));
if (batchSize <= 0)
{
break;
}
candidate = candidateSpan[index++];
powUsed++;
ReadOnlySpan<ulong> batch = candidateSpan.Slice(index, batchSize);
int startIndex = 0;
equalsOne = powStepper.ComputeNextIsUnity(candidate);
if (!equalsOne)
if (!powStepperInitialized)
{
continue;
ulong candidate = batch[0];
powUsed++;
bool equalsOne = powStepper.InitializeCpuIsUnity(candidate);
powStepperInitialized = true;
if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget))
{
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
startIndex = 1;
}
if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget))
for (int i = startIndex; i < batchSize && powUsed < powBudget; i++)
{
continue;
ulong candidate = batch[i];
powUsed++;
bool equalsOne = powStepper.ComputeNextIsUnity(candidate);
if (!equalsOne)
{
continue;
}
if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget))
{
continue;
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
result = candidate;
return true;
index += batchSize;
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
return false;
}
finally
@@ -572,7 +666,7 @@ internal static partial class PrimeOrderCalculator
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
{
int length = factors.Count;
// The phi factorization on the scanner path always yields at least one entry for phi >= 2,
@@ -582,13 +676,15 @@ internal static partial class PrimeOrderCalculator
// return true;
// }
ReadOnlySpan<ulong> factorSpan = new(factors.Factors, 0, length);
Span<ulong> stackBuffer = stackalloc ulong[length];
if (length <= 8)
{
return EvaluateSpecialMaxCandidatesCpu(factors, phi);
return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData);
}
ReadOnlySpan<ulong> factorSpan = new(factors.Factors, 0, length);
return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, divisorData);
return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData);
}
}
PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs
+140 -199
View File
@@ -1,90 +1,64 @@
using System.Buffers;
using System.Diagnostics;
using System.Numerics;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static void SortCandidatesWide(UInt128 prime, in UInt128? previousOrder, List<UInt128> candidates)
{
UInt128 previous = previousOrder ?? UInt128.Zero;
int previousGroup = previousOrder.HasValue ? GetGroupWide(previousOrder.Value, prime) : 1;
bool hasPrevious = previousOrder.HasValue;
int previousGroup = hasPrevious ? GetGroupWide(previous, prime) : 1;
candidates.Sort((left, right) => CompareCandidatesWide(prime, left, right, previous, hasPrevious, previousGroup));
candidates.Sort((x, y) =>
{
CandidateKey128 keyX = BuildKeyWide(x, prime, previous, previousGroup, hasPrevious);
CandidateKey128 keyY = BuildKeyWide(y, prime, previous, previousGroup, hasPrevious);
int primary = keyX.Primary.CompareTo(keyY.Primary);
if (primary != 0)
{
return primary;
}
int secondary = CompareComponents(keyX.SecondaryDescending, keyX.Secondary, keyY.SecondaryDescending, keyY.Secondary);
if (secondary != 0)
{
return secondary;
}
return CompareComponents(keyX.TertiaryDescending, keyX.Tertiary, keyY.TertiaryDescending, keyY.Tertiary);
});
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static int CompareCandidatesWide(UInt128 prime, UInt128 left, UInt128 right, UInt128 previous, bool hasPrevious, int previousGroup)
{
int leftGroup = GetGroupWide(left, prime);
int rightGroup = GetGroupWide(right, prime);
if (leftGroup == 0)
{
return rightGroup == 0 ? 0 : 1;
}
private static CandidateKey128 BuildKeyWide(
in UInt128 value,
in UInt128 prime,
in UInt128 previous,
int previousGroup,
bool hasPrevious)
{
int group = GetGroupWide(value, prime);
if (group == 0)
{
return new CandidateKey128(int.MaxValue, false, UInt128.Zero, false, UInt128.Zero);
}
if (rightGroup == 0)
{
return -1;
}
bool isGe = !hasPrevious || value >= previous;
int primary = ComputePrimary(group, isGe, previousGroup);
bool leftIsGe = !hasPrevious || left >= previous;
bool rightIsGe = !hasPrevious || right >= previous;
int leftPrimary = ComputePrimary(leftGroup, leftIsGe, previousGroup);
int rightPrimary = ComputePrimary(rightGroup, rightIsGe, previousGroup);
if (group == 3)
{
return new CandidateKey128(primary, true, value, true, value);
}
int compare = leftPrimary.CompareTo(rightPrimary);
if (compare != 0)
{
return compare;
}
UInt128 distance = hasPrevious ? (value > previous ? value - previous : previous - value) : value;
return new CandidateKey128(primary, false, distance, false, value);
}
bool leftDescending = leftGroup == 3;
bool rightDescending = rightGroup == 3;
UInt128 leftSecondary;
UInt128 leftTertiary;
if (leftDescending)
{
leftSecondary = left;
leftTertiary = left;
}
else
{
UInt128 leftReference = hasPrevious ? (left > previous ? left - previous : previous - left) : left;
leftSecondary = leftReference;
leftTertiary = left;
}
UInt128 rightSecondary;
UInt128 rightTertiary;
if (rightDescending)
{
rightSecondary = right;
rightTertiary = right;
}
else
{
UInt128 rightReference = hasPrevious ? (right > previous ? right - previous : previous - right) : right;
rightSecondary = rightReference;
rightTertiary = right;
}
compare = CompareComponents(leftDescending, leftSecondary, rightDescending, rightSecondary);
if (compare != 0)
{
return compare;
}
return CompareComponents(leftDescending, leftTertiary, rightDescending, rightTertiary);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY)
private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY)
{
if (descendingX == descendingY)
{
@@ -109,23 +83,22 @@ internal static partial class PrimeOrderCalculator
return valueX > valueY ? 1 : -1;
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static int GetGroupWide(in UInt128 value, in UInt128 prime)
{
UInt128 threshold = prime >> 3;
if (value <= threshold)
UInt128 threshold1 = prime >> 3;
if (value <= threshold1)
{
return 1;
}
threshold = prime >> 2;
if (value <= threshold)
UInt128 threshold2 = prime >> 2;
if (value <= threshold2)
{
return 2;
}
threshold = (UInt128)(((BigInteger)prime * 3) >> 3);
if (value <= threshold)
UInt128 threshold3 = (UInt128)(((BigInteger)prime * 3) >> 3);
if (value <= threshold3)
{
return 3;
}
@@ -133,10 +106,8 @@ internal static partial class PrimeOrderCalculator
return 0;
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List<UInt128> candidates, int limit)
{
// TODO: Is this condition ever met on EvenPerfectBitScanner execution paths for this method? Can we remove it?
if (count == 0)
{
return;
@@ -147,7 +118,6 @@ internal static partial class PrimeOrderCalculator
BuildCandidatesRecursiveWide(order, buffer, 0, UInt128.One, candidates, limit);
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static void BuildCandidatesRecursiveWide(
in UInt128 order,
in ReadOnlySpan<FactorEntry128> factors,
@@ -215,14 +185,33 @@ internal static partial class PrimeOrderCalculator
return PartialFactorResult128.Empty;
}
PartialFactorResult128 result = PartialFactorResult128.Rent();
Dictionary<UInt128, int> counts = result.FactorCounts;
counts.Clear();
Dictionary<UInt128, int> counts = new(capacity: 8);
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
UInt128 remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
UInt128 remaining;
List<UInt128> pending = result.Pending;
pending.Clear();
// if (IsGpuPow2Allowed && value <= ulong.MaxValue)
// {
// Dictionary<ulong, int> narrowCounts = new(capacity: 8);
// if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining))
// {
// foreach (KeyValuePair<ulong, int> entry in narrowCounts)
// {
// counts[(UInt128)entry.Key] = entry.Value;
// }
// remaining = (UInt128)narrowRemaining;
// }
// else
// {
// remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
// }
// }
// else
// {
remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
// }
List<UInt128> pending = new();
if (remaining > UInt128.One)
{
pending.Add(remaining);
@@ -231,12 +220,11 @@ internal static partial class PrimeOrderCalculator
if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
{
long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
FixedCapacityStack<UInt128> stack = result.CompositeStack;
stack.Clear();
Stack<UInt128> stack = new();
stack.Push(remaining);
pending.Clear();
long timestamp; // reused for deadline checks.
long timestamp = 0L; // reused for deadline checks.
while (stack.Count > 0)
{
UInt128 composite = stack.Pop();
@@ -288,13 +276,33 @@ internal static partial class PrimeOrderCalculator
if (counts.Count == 0 && cofactor == value)
{
result.Cofactor = value;
// Everything else stays at default values initiated in .Rent.
return result;
return new PartialFactorResult128(null, value, false, 0);
}
result.InitializeFromCounts(counts, cofactor, cofactor == UInt128.One);
return result;
ArrayPool<FactorEntry128>? pool = null;
FactorEntry128[]? rented = counts.Count > 0 ? (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null;
int index = 0;
if (rented is not null)
{
foreach (KeyValuePair<UInt128, int> entry in counts)
{
rented[index] = new FactorEntry128(entry.Key, entry.Value);
index++;
}
}
FactorEntry128[]? resultArray = null;
if (pool is not null)
{
Span<FactorEntry128> span = rented.AsSpan(0, index);
span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
resultArray = new FactorEntry128[index];
span.CopyTo(resultArray);
pool.Return(rented!, clearArray: false);
}
bool fullyFactored = cofactor == UInt128.One;
return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index);
}
private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary<UInt128, int> counts)
@@ -349,7 +357,7 @@ internal static partial class PrimeOrderCalculator
UInt128 c = UInt128.One;
UInt128 x = (UInt128)2UL;
UInt128 y = x;
long timestamp; // reused for deadline checks.
long timestamp = 0L; // reused for deadline checks.
while (true)
{
@@ -360,7 +368,8 @@ internal static partial class PrimeOrderCalculator
}
x = AdvancePolynomialWide(x, c, n);
y = AdvancePolynomialTwiceWide(y, c, n);
y = AdvancePolynomialWide(y, c, n);
y = AdvancePolynomialWide(y, c, n);
UInt128 diff = x > y ? x - y : y - x;
UInt128 d = BinaryGcdWide(diff, n);
@@ -387,14 +396,6 @@ internal static partial class PrimeOrderCalculator
return (UInt128)value;
}
private static UInt128 AdvancePolynomialTwiceWide(in UInt128 x, in UInt128 c, in UInt128 modulus)
{
BigInteger value = (BigInteger)x;
value = (value * value + (BigInteger)c) % (BigInteger)modulus;
value = (value * value + (BigInteger)c) % (BigInteger)modulus;
return (UInt128)value;
}
private static UInt128 BinaryGcdWide(UInt128 a, UInt128 b)
{
if (a == UInt128.Zero)
@@ -478,7 +479,8 @@ internal static partial class PrimeOrderCalculator
while (d == UInt128.One)
{
x = AdvancePolynomialWide(x, c, n);
y = AdvancePolynomialTwiceWide(y, c, n);
y = AdvancePolynomialWide(y, c, n);
y = AdvancePolynomialWide(y, c, n);
UInt128 diff = x > y ? x - y : y - x;
d = BinaryGcdWide(diff, n);
@@ -507,13 +509,12 @@ internal static partial class PrimeOrderCalculator
}
uint[] smallPrimes = PrimesGenerator.SmallPrimes;
ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
int length = smallPrimes.Length;
for (int i = 0; i < length; i++)
{
uint prime = smallPrimes[i];
UInt128 primeValue = prime;
UInt128 primeSquare = (UInt128)smallPrimesPow2[i];
UInt128 primeSquare = primeValue * primeValue;
if (primeSquare > value)
{
break;
@@ -587,104 +588,44 @@ internal static partial class PrimeOrderCalculator
}
}
private sealed class PartialFactorResult128
private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary)
{
[ThreadStatic]
private static PartialFactorResult128? s_poolHead;
public int Primary { get; } = primary;
public static readonly PartialFactorResult128 Empty = new()
{
Cofactor = UInt128.One,
FullyFactored = true,
Count = 0
};
public bool SecondaryDescending { get; } = secondaryDescending;
private const int ExponentHardLimit = 256;
public UInt128 Secondary { get; } = secondary;
internal readonly FactorEntry128[] TempFactors = new FactorEntry128[ExponentHardLimit];
internal readonly List<UInt128> FactorCandidatesList = new(ExponentHardLimit);
internal readonly Dictionary<UInt128, int> FactorCounts = new(capacity: ExponentHardLimit);
internal readonly List<UInt128> Pending = new(capacity: ExponentHardLimit);
internal readonly FixedCapacityStack<UInt128> CompositeStack = new(ExponentHardLimit);
public bool TertiaryDescending { get; } = tertiaryDescending;
private PartialFactorResult128? _next;
public UInt128 Cofactor;
public int Count;
public readonly FactorEntry128[] Factors = new FactorEntry128[ExponentHardLimit];
public bool FullyFactored;
public bool HasFactors;
private PartialFactorResult128()
{
}
public static PartialFactorResult128 Rent()
{
if (s_poolHead is { } instance)
{
s_poolHead = instance._next;
instance._next = null;
instance.Cofactor = UInt128.One;
instance.Count = 0;
instance.FullyFactored = false;
instance.HasFactors = false;
return instance;
}
return new PartialFactorResult128();
}
public void WithAdditionalPrime(UInt128 prime)
{
int index = Count;
FactorEntry128[] buffer = Factors;
buffer[index] = new FactorEntry128(prime, 1);
index++;
Count = index;
buffer.AsSpan(0, index).Sort(static (a, b) => a.Value.CompareTo(b.Value));
Cofactor = UInt128.One;
FullyFactored = true;
HasFactors = true;
}
public void Dispose()
{
if (this == Empty)
{
return;
}
_next = s_poolHead;
s_poolHead = this;
}
internal void InitializeFromCounts(Dictionary<UInt128, int> counts, UInt128 cofactor, bool fullyFactored)
{
int count = counts.Count;
Count = count;
Cofactor = cofactor;
FullyFactored = fullyFactored;
bool hasFactors = count != 0;
HasFactors = hasFactors;
if (!hasFactors)
{
return;
}
int index = 0;
foreach (KeyValuePair<UInt128, int> entry in counts)
{
Factors[index++] = new FactorEntry128(entry.Key, entry.Value);
}
// We benefit from static (non-virtual) call to Sort. When JIT finds a method call after a constructor, it optimizes it to a true static call.
new Span<FactorEntry128>(Factors, 0, count).Sort(static (a, b) => a.Value.CompareTo(b.Value));
}
public UInt128 Tertiary { get; } = tertiary;
}
private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count)
{
public FactorEntry128[]? Factors { get; } = factors;
public UInt128 Cofactor { get; } = cofactor;
public bool FullyFactored { get; } = fullyFactored;
public int Count { get; } = count;
public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
public PartialFactorResult128 WithAdditionalPrime(UInt128 prime)
{
if (Factors is null)
{
FactorEntry128[] local = new FactorEntry128[1];
local[0] = new FactorEntry128(prime, 1);
return new PartialFactorResult128(local, UInt128.One, true, 1);
}
FactorEntry128[] extended = new FactorEntry128[Count + 1];
Array.Copy(Factors, extended, Count);
extended[Count] = new FactorEntry128(prime, 1);
return new PartialFactorResult128(extended, UInt128.One, true, Count + 1);
}
}
}
PerfectNumbers.Core/ThreadStaticPools.cs
+5 -5
View File
@@ -263,12 +263,12 @@ namespace PerfectNumbers.Core
public static void Return((ulong, ulong)[] array) => _ulongUlongTupleArrayPool!.Return(array, clearArray: false);
[ThreadStatic]
private static FixedCapacityStack<List<PartialFactorPendingEntry>>? _primeOrderPendingEntryListPool;
private static FixedCapacityStack<List<PrimeOrderCalculator.PendingEntry>>? _primeOrderPendingEntryListPool;
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
internal static List<PartialFactorPendingEntry> RentPrimeOrderPendingEntryList(int capacityHint)
internal static List<PrimeOrderCalculator.PendingEntry> RentPrimeOrderPendingEntryList(int capacityHint)
{
FixedCapacityStack<List<PartialFactorPendingEntry>>? pool = _primeOrderPendingEntryListPool;
FixedCapacityStack<List<PrimeOrderCalculator.PendingEntry>>? pool = _primeOrderPendingEntryListPool;
if (pool is not null)
{
if (pool.Pop() is { } list)
@@ -286,11 +286,11 @@ namespace PerfectNumbers.Core
_primeOrderPendingEntryListPool = new(PerfectNumberConstants.DefaultPoolCapacity);
}
return new List<PartialFactorPendingEntry>(capacityHint);
return new List<PrimeOrderCalculator.PendingEntry>(capacityHint);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
internal static void ReturnPrimeOrderPendingEntryList(List<PartialFactorPendingEntry> list) => _primeOrderPendingEntryListPool!.Push(list);
internal static void ReturnPrimeOrderPendingEntryList(List<PrimeOrderCalculator.PendingEntry> list) => _primeOrderPendingEntryListPool!.Push(list);
[ThreadStatic]
private static FixedCapacityStack<FixedCapacityStack<ulong>>? _ulongStackPool;