Batch CPU divisor scanning across reusable buffers
This commit is contained in:
Marek Lemańczyk committed
1 parent
e5486af49d
commit
681e9a4a5e
1 file changed
+308
-84
@@ -1,6 +1,6 @@
|
||||
using System;
|
||||
using System.Buffers;
|
||||
using System.Collections.Concurrent;
|
||||
using System.Collections.Generic;
|
||||
using System.Numerics;
|
||||
using PerfectNumbers.Core;
|
||||
|
||||
@@ -200,12 +200,25 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
return false;
|
||||
}
|
||||
|
||||
Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
|
||||
int batchCapacity = Math.Max(1, _batchSize);
|
||||
DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
|
||||
int preferredCycleBatch = cycleCache.PreferredBatchSize;
|
||||
if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity)
|
||||
{
|
||||
batchCapacity = preferredCycleBatch;
|
||||
}
|
||||
|
||||
ulong cachedModulus = 0UL;
|
||||
MontgomeryDivisorData cachedDivisorData = default;
|
||||
bool hasCachedDivisorData = false;
|
||||
ulong[] candidateBuffer = ArrayPool<ulong>.Shared.Rent(batchCapacity);
|
||||
ulong[] cycleBuffer = ArrayPool<ulong>.Shared.Rent(batchCapacity);
|
||||
MontgomeryDivisorData[] divisorDataBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity);
|
||||
byte[] remainderBuffer = ArrayPool<byte>.Shared.Rent(batchCapacity * 6);
|
||||
|
||||
Span<byte> remainder10Span = remainderBuffer.AsSpan(0, batchCapacity);
|
||||
Span<byte> remainder8Span = remainderBuffer.AsSpan(batchCapacity, batchCapacity);
|
||||
Span<byte> remainder5Span = remainderBuffer.AsSpan(batchCapacity * 2, batchCapacity);
|
||||
Span<byte> remainder3Span = remainderBuffer.AsSpan(batchCapacity * 3, batchCapacity);
|
||||
Span<byte> remainder7Span = remainderBuffer.AsSpan(batchCapacity * 4, batchCapacity);
|
||||
Span<byte> remainder11Span = remainderBuffer.AsSpan(batchCapacity * 5, batchCapacity);
|
||||
|
||||
byte step10 = (byte)(step % 10UL);
|
||||
byte step8 = (byte)(step % 8UL);
|
||||
@@ -221,83 +234,151 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
byte remainder7 = (byte)(divisor % 7UL);
|
||||
byte remainder11 = (byte)(divisor % 11UL);
|
||||
|
||||
// Keep the divisibility filters aligned with the divisor-cycle generator so the
|
||||
// CPU path never requests cycles that were skipped during cache creation.
|
||||
bool lastIsSeven = (prime & 3UL) == 3UL;
|
||||
|
||||
while (divisor <= limit)
|
||||
try
|
||||
{
|
||||
if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired())
|
||||
while (divisor <= limit)
|
||||
{
|
||||
processedAll = false;
|
||||
return false;
|
||||
}
|
||||
Span<ulong> candidateSpan = candidateBuffer.AsSpan(0, batchCapacity);
|
||||
int chunkCount = 0;
|
||||
UInt128 localDivisor = divisor;
|
||||
|
||||
iterationCounter++;
|
||||
ulong candidate = (ulong)divisor;
|
||||
processedCount++;
|
||||
lastProcessed = candidate;
|
||||
byte localRemainder10 = remainder10;
|
||||
byte localRemainder8 = remainder8;
|
||||
byte localRemainder5 = remainder5;
|
||||
byte localRemainder3 = remainder3;
|
||||
byte localRemainder7 = remainder7;
|
||||
byte localRemainder11 = remainder11;
|
||||
|
||||
bool admissible = lastIsSeven
|
||||
? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9)
|
||||
: (remainder10 == 1 || remainder10 == 3 || remainder10 == 9);
|
||||
|
||||
if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0)
|
||||
{
|
||||
if (!hasCachedDivisorData || cachedModulus != candidate)
|
||||
while (chunkCount < batchCapacity && localDivisor <= limit)
|
||||
{
|
||||
cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate);
|
||||
cachedModulus = candidate;
|
||||
hasCachedDivisorData = true;
|
||||
}
|
||||
|
||||
MontgomeryDivisorData divisorData = cachedDivisorData;
|
||||
ulong divisorCycle;
|
||||
|
||||
if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
|
||||
{
|
||||
divisorCycle = cycleCache.GetCycleLength(candidate);
|
||||
}
|
||||
else
|
||||
{
|
||||
factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
|
||||
// EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot.
|
||||
// Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator
|
||||
// can rely on an odd divisor and a prime exponent when deriving the cycle length.
|
||||
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
|
||||
if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired())
|
||||
{
|
||||
divisorCycle = cycleCache.GetCycleLength(candidate);
|
||||
processedAll = false;
|
||||
return false;
|
||||
}
|
||||
|
||||
ulong candidate = (ulong)localDivisor;
|
||||
candidateSpan[chunkCount] = candidate;
|
||||
remainder10Span[chunkCount] = localRemainder10;
|
||||
remainder8Span[chunkCount] = localRemainder8;
|
||||
remainder5Span[chunkCount] = localRemainder5;
|
||||
remainder3Span[chunkCount] = localRemainder3;
|
||||
remainder7Span[chunkCount] = localRemainder7;
|
||||
remainder11Span[chunkCount] = localRemainder11;
|
||||
|
||||
processedCount++;
|
||||
lastProcessed = candidate;
|
||||
iterationCounter++;
|
||||
chunkCount++;
|
||||
|
||||
localDivisor += step;
|
||||
localRemainder10 = AddMod(localRemainder10, step10, (byte)10);
|
||||
localRemainder8 = AddMod(localRemainder8, step8, (byte)8);
|
||||
localRemainder5 = AddMod(localRemainder5, step5, (byte)5);
|
||||
localRemainder3 = AddMod(localRemainder3, step3, (byte)3);
|
||||
localRemainder7 = AddMod(localRemainder7, step7, (byte)7);
|
||||
localRemainder11 = AddMod(localRemainder11, step11, (byte)11);
|
||||
}
|
||||
|
||||
if (chunkCount == 0)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
divisor = localDivisor;
|
||||
remainder10 = localRemainder10;
|
||||
remainder8 = localRemainder8;
|
||||
remainder5 = localRemainder5;
|
||||
remainder3 = localRemainder3;
|
||||
remainder7 = localRemainder7;
|
||||
remainder11 = localRemainder11;
|
||||
|
||||
// Reuse the candidate buffer in-place so the filtered results stay contiguous without additional arrays.
|
||||
Span<ulong> chunkCandidates = candidateSpan[..chunkCount];
|
||||
int filteredCount = 0;
|
||||
|
||||
for (int i = 0; i < chunkCount; i++)
|
||||
{
|
||||
bool admissible = lastIsSeven
|
||||
? (remainder10Span[i] == 3 || remainder10Span[i] == 7 || remainder10Span[i] == 9)
|
||||
: (remainder10Span[i] == 1 || remainder10Span[i] == 3 || remainder10Span[i] == 9);
|
||||
|
||||
if (!admissible)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
if ((remainder8Span[i] != 1 && remainder8Span[i] != 7)
|
||||
|| remainder3Span[i] == 0
|
||||
|| remainder5Span[i] == 0
|
||||
|| remainder7Span[i] == 0
|
||||
|| remainder11Span[i] == 0)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
// Store the accepted candidate at the next free slot so the front of the span becomes the filtered set.
|
||||
chunkCandidates[filteredCount++] = chunkCandidates[i];
|
||||
}
|
||||
|
||||
if (filteredCount == 0)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
Span<ulong> filteredSpan = chunkCandidates[..filteredCount];
|
||||
Span<ulong> cycleSpan = cycleBuffer.AsSpan(0, filteredCount);
|
||||
cycleCache.GetCycleLengths(filteredSpan, cycleSpan);
|
||||
|
||||
int admissibleCount = 0;
|
||||
Span<MontgomeryDivisorData> divisorDataSpan = divisorDataBuffer.AsSpan(0, filteredCount);
|
||||
|
||||
for (int i = 0; i < filteredCount; i++)
|
||||
{
|
||||
if (cycleSpan[i] != prime)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
ulong candidate = filteredSpan[i];
|
||||
// Reuse the filtered span's leading segment to accumulate the admissible divisors for Montgomery checks.
|
||||
filteredSpan[admissibleCount] = candidate;
|
||||
divisorDataSpan[admissibleCount] = MontgomeryDivisorData.FromModulus(candidate);
|
||||
admissibleCount++;
|
||||
}
|
||||
|
||||
if (admissibleCount == 0)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
Span<ulong> admissibleSpan = filteredSpan[..admissibleCount];
|
||||
Span<MontgomeryDivisorData> admissibleDataSpan = divisorDataSpan[..admissibleCount];
|
||||
|
||||
for (int i = 0; i < admissibleCount; i++)
|
||||
{
|
||||
if (CheckDivisor(prime, prime, in admissibleDataSpan[i]) != 0)
|
||||
{
|
||||
lastProcessed = admissibleSpan[i];
|
||||
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;
|
||||
remainder10 = AddMod(remainder10, step10, (byte)10);
|
||||
remainder8 = AddMod(remainder8, step8, (byte)8);
|
||||
remainder5 = AddMod(remainder5, step5, (byte)5);
|
||||
remainder3 = AddMod(remainder3, step3, (byte)3);
|
||||
remainder7 = AddMod(remainder7, step7, (byte)7);
|
||||
remainder11 = AddMod(remainder11, step11, (byte)11);
|
||||
}
|
||||
finally
|
||||
{
|
||||
ArrayPool<ulong>.Shared.Return(candidateBuffer, clearArray: false);
|
||||
ArrayPool<ulong>.Shared.Return(cycleBuffer, clearArray: false);
|
||||
ArrayPool<MontgomeryDivisorData>.Shared.Return(divisorDataBuffer, clearArray: false);
|
||||
ArrayPool<byte>.Shared.Return(remainderBuffer, clearArray: false);
|
||||
}
|
||||
|
||||
processedAll = divisor > limit;
|
||||
return false;
|
||||
}
|
||||
|
||||
private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
|
||||
{
|
||||
ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
|
||||
@@ -372,6 +453,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
private bool _hasCachedDivisorData;
|
||||
private ulong _cachedDivisor;
|
||||
private MontgomeryDivisorData _cachedDivisorData;
|
||||
private ulong[]? _primeDeltas;
|
||||
private ulong[]? _cycleRemainders;
|
||||
private ulong[]? _montgomeryDeltas;
|
||||
private int _bufferCapacity;
|
||||
|
||||
internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner)
|
||||
{
|
||||
@@ -412,6 +497,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
cachedData = _cachedDivisorData;
|
||||
}
|
||||
|
||||
ulong modulus = cachedData.Modulus;
|
||||
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
|
||||
{
|
||||
hits.Clear();
|
||||
return;
|
||||
}
|
||||
|
||||
if (divisorCycle == 0UL)
|
||||
{
|
||||
divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
|
||||
@@ -422,35 +514,167 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
}
|
||||
}
|
||||
|
||||
// Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
|
||||
// They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
|
||||
var exponentStepper = new ExponentRemainderStepper(cachedData);
|
||||
if (!exponentStepper.IsValidModulus)
|
||||
{
|
||||
hits.Clear();
|
||||
return;
|
||||
}
|
||||
EnsureCapacity(length);
|
||||
|
||||
var cycleStepper = new CycleRemainderStepper(divisorCycle);
|
||||
Span<ulong> deltaSpan = _primeDeltas!.AsSpan(0, length);
|
||||
ComputePrimeDeltas(primes, deltaSpan);
|
||||
|
||||
ulong remainder = cycleStepper.Initialize(primes[0]);
|
||||
hits[0] = remainder == 0UL
|
||||
? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0)
|
||||
: (byte)0;
|
||||
Span<ulong> remainderSpan = _cycleRemainders!.AsSpan(0, length);
|
||||
ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
|
||||
|
||||
Span<ulong> montgomerySpan = _montgomeryDeltas!.AsSpan(0, length);
|
||||
ComputeMontgomeryDeltas(cachedData, primes, deltaSpan, montgomerySpan);
|
||||
|
||||
ulong nPrime = cachedData.NPrime;
|
||||
ulong montgomeryOne = cachedData.MontgomeryOne;
|
||||
ulong currentMontgomery = montgomerySpan[0];
|
||||
|
||||
hits[0] = remainderSpan[0] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0;
|
||||
|
||||
for (int i = 1; i < length; i++)
|
||||
{
|
||||
remainder = cycleStepper.ComputeNext(primes[i]);
|
||||
if (remainder != 0UL)
|
||||
currentMontgomery = currentMontgomery.MontgomeryMultiply(montgomerySpan[i], modulus, nPrime);
|
||||
hits[i] = remainderSpan[i] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0;
|
||||
}
|
||||
}
|
||||
|
||||
private void EnsureCapacity(int requiredLength)
|
||||
{
|
||||
if (requiredLength <= 0)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
int newCapacity = _bufferCapacity;
|
||||
if (newCapacity == 0)
|
||||
{
|
||||
newCapacity = 1;
|
||||
}
|
||||
|
||||
while (newCapacity < requiredLength)
|
||||
{
|
||||
newCapacity <<= 1;
|
||||
}
|
||||
|
||||
ReturnBuffers();
|
||||
|
||||
_primeDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
|
||||
_cycleRemainders = ArrayPool<ulong>.Shared.Rent(newCapacity);
|
||||
_montgomeryDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
|
||||
_bufferCapacity = newCapacity;
|
||||
}
|
||||
|
||||
private void ReturnBuffers()
|
||||
{
|
||||
if (_primeDeltas is not null)
|
||||
{
|
||||
ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
|
||||
_primeDeltas = null;
|
||||
}
|
||||
|
||||
if (_cycleRemainders is not null)
|
||||
{
|
||||
ArrayPool<ulong>.Shared.Return(_cycleRemainders, clearArray: false);
|
||||
_cycleRemainders = null;
|
||||
}
|
||||
|
||||
if (_montgomeryDeltas is not null)
|
||||
{
|
||||
ArrayPool<ulong>.Shared.Return(_montgomeryDeltas, clearArray: false);
|
||||
_montgomeryDeltas = null;
|
||||
}
|
||||
|
||||
_bufferCapacity = 0;
|
||||
}
|
||||
|
||||
private static void ComputePrimeDeltas(ReadOnlySpan<ulong> primes, Span<ulong> deltas)
|
||||
{
|
||||
if (primes.Length == 0)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
ulong previous = primes[0];
|
||||
deltas[0] = previous;
|
||||
|
||||
for (int i = 1; i < primes.Length; i++)
|
||||
{
|
||||
ulong current = primes[i];
|
||||
if (current <= previous)
|
||||
{
|
||||
hits[i] = 0;
|
||||
continue;
|
||||
throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing.");
|
||||
}
|
||||
|
||||
hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0;
|
||||
deltas[i] = current - previous;
|
||||
previous = current;
|
||||
}
|
||||
}
|
||||
|
||||
private static void ComputeCycleRemainders(ulong cycleLength, ReadOnlySpan<ulong> primes, ReadOnlySpan<ulong> deltas, Span<ulong> remainders)
|
||||
{
|
||||
if (cycleLength == 0UL || primes.Length == 0)
|
||||
{
|
||||
remainders.Clear();
|
||||
return;
|
||||
}
|
||||
|
||||
ulong remainder = primes[0] % cycleLength;
|
||||
remainders[0] = remainder;
|
||||
|
||||
for (int i = 1; i < primes.Length; i++)
|
||||
{
|
||||
remainder = AddMod(remainder, deltas[i], cycleLength);
|
||||
remainders[i] = remainder;
|
||||
}
|
||||
}
|
||||
|
||||
private static void ComputeMontgomeryDeltas(in MontgomeryDivisorData divisorData, ReadOnlySpan<ulong> primes, ReadOnlySpan<ulong> deltas, Span<ulong> montgomeryValues)
|
||||
{
|
||||
if (primes.Length == 0)
|
||||
{
|
||||
montgomeryValues.Clear();
|
||||
return;
|
||||
}
|
||||
|
||||
montgomeryValues[0] = primes[0].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true);
|
||||
|
||||
for (int i = 1; i < primes.Length; i++)
|
||||
{
|
||||
montgomeryValues[i] = deltas[i].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true);
|
||||
}
|
||||
}
|
||||
|
||||
private static ulong AddMod(ulong value, ulong delta, ulong modulus)
|
||||
{
|
||||
if (modulus == 0UL)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
|
||||
if (ulong.MaxValue - value < delta)
|
||||
{
|
||||
UInt128 extended = (UInt128)value + delta;
|
||||
return (ulong)(extended % modulus);
|
||||
}
|
||||
|
||||
ulong sum = value + delta;
|
||||
if (sum >= modulus)
|
||||
{
|
||||
sum -= modulus;
|
||||
if (sum >= modulus)
|
||||
{
|
||||
sum %= modulus;
|
||||
}
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
public void Dispose()
|
||||
{
|
||||
if (_disposed)
|
||||
|
||||
Reference in new issue
Block a user