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

Merge pull request #218 from mlemanczyk/codex/analyze-mersenne=bydivisor-path-for-tracking-remainder-bbul2g

Browse Source 
Reuse stack spans across divisor processing loops
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-10-21 20:54:45 +02:00
2 parents 31a0b8d706 a226ae10cb
commit 9bcf1bd9a1
4 files changed
+530 -147
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
+326 -54
View File
@@ -1,3 +1,4 @@
using System;
using System.Buffers;
using System.Diagnostics;
using System.Numerics;
@@ -110,27 +111,70 @@ internal static partial class PrimeOrderCalculator
private static bool TrySpecialMaxCpu(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,
// preventing the zero-length case from occurring outside synthetic tests.
// if (length == 0)
// {
// return true;
// }
ReadOnlySpan<FactorEntry> factorSpan = new(factors.Factors, 0, length);
for (int i = 0; i < length; i++)
if (length <= 32)
{
ulong factor = factorSpan[i].Value;
// Don't uncomment this. If there is division by zero, that means something wrong.
// if (factor == 0UL)
Span<ulong> stackBuffer = stackalloc ulong[length];
return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData);
}
ulong[] rented = ThreadStaticPools.UlongPool.Rent(length);
bool result = EvaluateSpecialMaxCandidates(rented.AsSpan(0, length), factorSpan, phi, prime, divisorData);
ThreadStaticPools.UlongPool.Return(rented, clearArray: false);
return result;
}
private static bool EvaluateSpecialMaxCandidates(Span<ulong> buffer, ReadOnlySpan<FactorEntry> factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData)
{
int actual = 0;
int factorCount = factors.Length;
for (int i = 0; i < factorCount; i++)
{
ulong factor = factors[i].Value;
// The partial factor pipeline never feeds zero or oversized factors while scanning candidate orders.
// if (factor == 0UL || factor > phi)
// {
// continue;
// }
// TODO: Can we optimize it thanks to knowing the factors, differences between factors, factors and / or phi
// being monotonically increasing etc. For example, if phi increases and factor increases, what does it tell us
// about the value of reduced? Or, maybe wee can use residue stepping to calculate the residue instead of calling
// Pow2equalsOneCpu? Maybe we don't need to calculate it from scratch?
ulong reduced = phi / factor;
if (Pow2EqualsOneCpu(reduced, prime, divisorData))
if (reduced == 0UL)
{
continue;
}
buffer[actual] = reduced;
actual++;
}
if (actual == 0)
{
return true;
}
Span<ulong> candidates = buffer.Slice(0, actual);
candidates.Sort();
ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
int candidateCount = candidates.Length;
for (int i = 0; i < candidateCount; i++)
{
if (stepper.ComputeNextIsUnity(candidates[i]))
{
ThreadStaticPools.ReturnExponentStepper(stepper);
return false;
}
}
ThreadStaticPools.ReturnExponentStepper(stepper);
return true;
}
@@ -160,8 +204,6 @@ internal static partial class PrimeOrderCalculator
Span<FactorEntry> buffer = tempArray;
factorSpan.CopyTo(buffer);
// bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor);
if (!factors.FullyFactored && factors.Cofactor > 1UL && factors.CofactorIsPrime)
{
buffer[length] = new FactorEntry(factors.Cofactor, 1, true);
@@ -170,32 +212,110 @@ internal static partial class PrimeOrderCalculator
buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value));
ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
const int StackExponentCapacity = 16;
const int ExponentHardLimit = 256;
Span<ulong> stackCandidates = stackalloc ulong[StackExponentCapacity];
Span<bool> stackEvaluations = stackalloc bool[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
bool[]? heapEvaluationArray = null;
for (int i = 0; i < length; i++)
{
ulong primeFactor = buffer[i].Value;
int exponent = buffer[i].Exponent;
for (int iteration = 0; iteration < exponent; iteration++)
// Factor exponents produced by Pollard-Rho and trial division are always positive in the scanner flow.
// if (exponent <= 0)
// {
// continue;
// }
if (exponent > ExponentHardLimit)
{
if ((order % primeFactor) != 0UL)
{
break;
}
ulong reduced = order / primeFactor;
if (Pow2EqualsOneCpu(reduced, prime, divisorData))
{
order = reduced;
continue;
}
break;
throw new InvalidOperationException($"Prime factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}.");
}
if (exponent <= StackExponentCapacity)
{
ProcessExponentLoweringPrime(stackCandidates.Slice(0, exponent), stackEvaluations.Slice(0, exponent), ref order, primeFactor, exponent, ref stepper);
continue;
}
heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit);
heapEvaluationArray ??= ThreadStaticPools.BoolPool.Rent(ExponentHardLimit);
ProcessExponentLoweringPrime(heapCandidateArray.AsSpan(0, exponent), heapEvaluationArray.AsSpan(0, exponent), ref order, primeFactor, exponent, ref stepper);
}
ThreadStaticPools.ReturnExponentStepper(stepper);
if (heapCandidateArray is not null)
{
ThreadStaticPools.UlongPool.Return(heapCandidateArray);
ThreadStaticPools.BoolPool.Return(heapEvaluationArray!);
}
pool.Return(tempArray, clearArray: false);
return order;
}
private static void ProcessExponentLoweringPrime(Span<ulong> candidateBuffer, Span<bool> evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepper stepper)
{
ulong working = order;
int actual = 0;
while (actual < exponent)
{
if (working < primeFactor)
{
break;
}
ulong reduced = working / primeFactor;
if (reduced == 0UL)
{
break;
}
candidateBuffer[actual] = reduced;
working = reduced;
actual++;
}
if (actual == 0)
{
return;
}
Span<ulong> candidates = candidateBuffer.Slice(0, actual);
Span<bool> evaluations = evaluationBuffer.Slice(0, actual);
for (int left = 0, right = actual - 1; left < right; left++, right--)
{
ulong tmp = candidates[left];
candidates[left] = candidates[right];
candidates[right] = tmp;
}
stepper.Reset();
for (int j = 0; j < actual; j++)
{
evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]);
}
for (int j = actual - 1; j >= 0; j--)
{
if (!evaluations[j])
{
break;
}
order = candidates[j];
}
}
private static bool TryConfirmOrderCpu(
ulong prime,
ulong order,
@@ -253,28 +373,101 @@ internal static partial class PrimeOrderCalculator
{
ReadOnlySpan<FactorEntry> span = factorization.Factors!;
int length = factorization.Count;
const int StackExponentCapacity = 32;
const int ExponentHardLimit = 256;
Span<ulong> stackBuffer = stackalloc ulong[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
bool violates = false;
for (int i = 0; i < length; i++)
{
ulong primeFactor = span[i].Value;
ulong reduced = order;
int exponent = span[i].Exponent;
for (int iteration = 0; iteration < exponent; iteration++)
// Factor exponents emitted by TryConfirmCandidateCpu stay positive for production workloads.
// if (exponent <= 0)
// {
// continue;
// }
if (exponent > ExponentHardLimit)
{
// TODO: Benchmark and / or implement DivRem like solution or residue stepper?
if ((reduced % primeFactor) != 0UL)
throw new InvalidOperationException($"Factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}.");
}
if (exponent <= StackExponentCapacity)
{
if (ValidateOrderForFactor(stackBuffer.Slice(0, exponent), primeFactor, exponent, order, ref stepper))
{
violates = true;
break;
}
reduced /= primeFactor;
if (Pow2EqualsOneCpu(reduced, prime, divisorData))
{
return false;
}
continue;
}
heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit);
if (ValidateOrderForFactor(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, order, ref stepper))
{
violates = true;
break;
}
}
return true;
ThreadStaticPools.ReturnExponentStepper(stepper);
if (heapCandidateArray is not null)
{
ThreadStaticPools.UlongPool.Return(heapCandidateArray);
}
return !violates;
}
private static bool ValidateOrderForFactor(Span<ulong> buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepper stepper)
{
ulong working = order;
int actual = 0;
while (actual < exponent)
{
if (working < primeFactor)
{
break;
}
ulong reduced = working / primeFactor;
if (reduced == 0UL)
{
break;
}
buffer[actual] = reduced;
working = reduced;
actual++;
}
if (actual == 0)
{
return false;
}
Span<ulong> candidates = buffer.Slice(0, actual);
candidates.Reverse();
stepper.Reset();
for (int j = 0; j < actual; j++)
{
if (stepper.ComputeNextIsUnity(candidates[j]))
{
return true;
}
}
return false;
}
private static bool TryHeuristicFinishCpu(
@@ -328,6 +521,7 @@ internal static partial class PrimeOrderCalculator
BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
if (candidates.Count == 0)
{
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
return false;
}
@@ -340,6 +534,7 @@ internal static partial class PrimeOrderCalculator
int candidateCount = candidates.Count;
bool allowGpuBatch = true;
Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
ExponentRemainderStepper powStepper = ThreadStaticPools.RentExponentStepper(divisorData);
// DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
int index = 0;
@@ -412,7 +607,7 @@ internal static partial class PrimeOrderCalculator
}
else
{
equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData);
equalsOne = powStepper.ComputeNextIsUnity(candidate);
}
if (!equalsOne)
{
@@ -429,7 +624,9 @@ internal static partial class PrimeOrderCalculator
pool.Return(gpuPool, clearArray: false);
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepper(powStepper);
result = candidate;
return true;
}
@@ -442,7 +639,9 @@ internal static partial class PrimeOrderCalculator
index += batchSize;
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepper(powStepper);
// DebugLog("No candidate confirmed");
return false;
}
@@ -635,8 +834,6 @@ internal static partial class PrimeOrderCalculator
return false;
}
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
// bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
if (!factorization.CofactorIsPrime)
{
return false;
@@ -649,33 +846,57 @@ internal static partial class PrimeOrderCalculator
ReadOnlySpan<FactorEntry> span = factorization.Factors;
int length = factorization.Count;
const int StackExponentCapacity = 32;
const int ExponentHardLimit = 256;
Span<ulong> stackBuffer = stackalloc ulong[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
bool violates = false;
for (int i = 0; i < length; i++)
{
ulong primeFactor = span[i].Value;
ulong reduced = candidate;
for (int iteration = 0; iteration < span[i].Exponent; iteration++)
int exponent = span[i].Exponent;
// Factorization entries for candidate confirmation are always positive here.
// if (exponent <= 0)
// {
// continue;
// }
if (exponent > ExponentHardLimit)
{
// TODO: Measure and implement DivRem like solution or residue stepper to reuse calculation results from previous iteration
if ((reduced % primeFactor) != 0UL)
throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}.");
}
if (exponent <= StackExponentCapacity)
{
if (CheckCandidateViolation(stackBuffer.Slice(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
{
violates = true;
break;
}
reduced /= primeFactor;
if (powUsed >= powBudget && powBudget > 0)
{
return false;
}
continue;
}
powUsed++;
if (Pow2EqualsOneCpu(reduced, prime, divisorData))
{
return false;
}
heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit);
if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
{
violates = true;
break;
}
}
return true;
ThreadStaticPools.ReturnExponentStepper(stepper);
if (heapCandidateArray is not null)
{
ThreadStaticPools.UlongPool.Return(heapCandidateArray);
}
return !violates;
}
finally
{
@@ -683,6 +904,55 @@ internal static partial class PrimeOrderCalculator
}
}
private static bool CheckCandidateViolation(Span<ulong> buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepper stepper)
{
ulong working = candidate;
int actual = 0;
while (actual < exponent)
{
if (working < primeFactor)
{
break;
}
ulong reduced = working / primeFactor;
if (reduced == 0UL)
{
break;
}
buffer[actual] = reduced;
working = reduced;
actual++;
}
if (actual == 0)
{
return false;
}
Span<ulong> candidates = buffer.Slice(0, actual);
candidates.Reverse();
stepper.Reset();
for (int j = 0; j < actual; j++)
{
if (powUsed >= powBudget && powBudget > 0)
{
return true;
}
powUsed++;
if (stepper.ComputeNextIsUnity(candidates[j]))
{
return true;
}
}
return false;
}
private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
{
if (IsGpuPow2Allowed)
@@ -763,6 +1033,7 @@ internal static partial class PrimeOrderCalculator
}
factorCount = copyIndex;
counts.Clear();
ThreadStaticPools.ReturnUlongIntDictionary(counts);
counts = null;
}
@@ -967,6 +1238,7 @@ internal static partial class PrimeOrderCalculator
if (counts is not null)
{
counts.Clear();
ThreadStaticPools.ReturnUlongIntDictionary(counts);
}
@@ -996,7 +1268,7 @@ internal static partial class PrimeOrderCalculator
// value will never <= 1 in production code
// if (value <= 1UL)
// {
// return true;
// return true;
// }
int capacity = Math.Min(primeTargets.Length, exponentTargets.Length);