Add heuristic prime sieves for cached divisor sequences
This commit is contained in:
Marek Lemańczyk committed
1 parent
4ee6464130
commit
15c229665d
10 files changed
+1597
-364
No files matched your search
@@ -12,16 +12,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
private ulong _lastStatusDivisor;
|
||||
private int _batchSize = 1_024;
|
||||
|
||||
[ThreadStatic]
|
||||
private static GpuUInt128WorkSet _divisorScanGpuWorkSet;
|
||||
|
||||
private struct GpuUInt128WorkSet
|
||||
{
|
||||
public GpuUInt128 Step;
|
||||
public GpuUInt128 Divisor;
|
||||
public GpuUInt128 Limit;
|
||||
}
|
||||
|
||||
public int BatchSize
|
||||
{
|
||||
get => _batchSize;
|
||||
@@ -113,151 +103,37 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
out bool processedAll,
|
||||
out ulong processedCount)
|
||||
{
|
||||
processedCount = 0UL;
|
||||
processedAll = false;
|
||||
processedCount = 0UL;
|
||||
|
||||
// The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs.
|
||||
// Keeping the guard commented out documents the reasoning for benchmarks and tests.
|
||||
// if (allowedMax < 3UL)
|
||||
// {
|
||||
// return false;
|
||||
// }
|
||||
var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
|
||||
|
||||
ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet;
|
||||
|
||||
ref GpuUInt128 step = ref workSet.Step;
|
||||
step.High = 0UL;
|
||||
step.Low = prime;
|
||||
step.ShiftLeft(1);
|
||||
|
||||
ref GpuUInt128 limit = ref workSet.Limit;
|
||||
limit.High = 0UL;
|
||||
limit.Low = allowedMax;
|
||||
|
||||
ref GpuUInt128 divisor = ref workSet.Divisor;
|
||||
divisor.High = step.High;
|
||||
divisor.Low = step.Low;
|
||||
divisor.Add(1UL);
|
||||
if (divisor.CompareTo(limit) > 0)
|
||||
while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
|
||||
{
|
||||
processedAll = true;
|
||||
divisor.High = 0UL;
|
||||
divisor.Low = 0UL;
|
||||
step.High = 0UL;
|
||||
step.Low = 0UL;
|
||||
limit.High = 0UL;
|
||||
limit.Low = 0UL;
|
||||
return false;
|
||||
}
|
||||
processedCount = enumerator.ProcessedCount;
|
||||
|
||||
// Intentionally recomputes factorizations without a per-thread cache.
|
||||
// The previous factor cache recorded virtually no hits and only slowed down the scan.
|
||||
DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
|
||||
PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate);
|
||||
ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength(
|
||||
prime,
|
||||
in preparation,
|
||||
out _,
|
||||
out _,
|
||||
out _);
|
||||
|
||||
ulong stepHigh = step.High;
|
||||
ulong stepLow = step.Low;
|
||||
byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
|
||||
byte step8 = (byte)(stepLow % 8UL);
|
||||
byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL);
|
||||
byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL);
|
||||
byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL);
|
||||
byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL);
|
||||
|
||||
ulong divisorHigh = divisor.High;
|
||||
ulong divisorLow = divisor.Low;
|
||||
byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL);
|
||||
byte remainder8 = (byte)(divisorLow % 8UL);
|
||||
byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL);
|
||||
byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL);
|
||||
byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
|
||||
byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 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;
|
||||
const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
|
||||
const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
|
||||
|
||||
while (divisor.CompareTo(limit) <= 0)
|
||||
{
|
||||
ulong candidate = divisor.Low;
|
||||
processedCount++;
|
||||
|
||||
ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
|
||||
bool admissible = ((decimalMask >> remainder10) & 1) != 0;
|
||||
|
||||
if (admissible
|
||||
&& (remainder8 == 1 || remainder8 == 7)
|
||||
&& remainder3 != 0
|
||||
&& remainder5 != 0
|
||||
&& remainder7 != 0
|
||||
&& remainder11 != 0)
|
||||
if (divisorCycle == prime)
|
||||
{
|
||||
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
|
||||
ulong divisorCycle;
|
||||
// Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs.
|
||||
// if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
|
||||
// {
|
||||
// divisorCycle = cycleCache.GetCycleLength(candidate);
|
||||
// }
|
||||
// else
|
||||
{
|
||||
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
|
||||
candidate,
|
||||
prime,
|
||||
divisorData,
|
||||
out ulong computedCycle,
|
||||
out bool primeOrderFailed) || computedCycle == 0UL)
|
||||
{
|
||||
// Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles
|
||||
// for the exponents scanned here, so skip the unused cache fallback and compute directly.
|
||||
divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
|
||||
candidate,
|
||||
divisorData,
|
||||
skipPrimeOrderHeuristic: primeOrderFailed);
|
||||
}
|
||||
else
|
||||
{
|
||||
divisorCycle = computedCycle;
|
||||
}
|
||||
}
|
||||
|
||||
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.
|
||||
divisor.High = 0UL;
|
||||
divisor.Low = 0UL;
|
||||
step.High = 0UL;
|
||||
step.Low = 0UL;
|
||||
limit.High = 0UL;
|
||||
limit.Low = 0UL;
|
||||
processedAll = true;
|
||||
return true;
|
||||
}
|
||||
|
||||
if (divisorCycle == 0UL)
|
||||
{
|
||||
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
|
||||
}
|
||||
processedAll = true;
|
||||
return true;
|
||||
}
|
||||
|
||||
divisor.Add(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);
|
||||
if (divisorCycle == 0UL)
|
||||
{
|
||||
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
|
||||
}
|
||||
}
|
||||
|
||||
processedAll = divisor.CompareTo(limit) > 0;
|
||||
divisor.High = 0UL;
|
||||
divisor.Low = 0UL;
|
||||
step.High = 0UL;
|
||||
step.Low = 0UL;
|
||||
limit.High = 0UL;
|
||||
limit.Low = 0UL;
|
||||
processedCount = enumerator.ProcessedCount;
|
||||
processedAll = enumerator.Exhausted;
|
||||
return false;
|
||||
}
|
||||
|
||||
@@ -267,17 +143,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
||||
return residue == 1UL ? (byte)1 : (byte)0;
|
||||
}
|
||||
|
||||
private static byte AddMod(byte value, byte delta, byte modulus)
|
||||
{
|
||||
int sum = value + delta;
|
||||
if (sum >= modulus)
|
||||
{
|
||||
sum -= modulus;
|
||||
}
|
||||
|
||||
return (byte)sum;
|
||||
}
|
||||
|
||||
private void UpdateStatusUnsafe(ulong processedCount)
|
||||
{
|
||||
if (processedCount == 0UL)
|
||||
|
||||
@@ -41,6 +41,17 @@ internal static class PrimeTesterKernels
|
||||
results[index] = composite ? (byte)0 : (byte)1;
|
||||
}
|
||||
|
||||
public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView<ulong> divisors, ulong n, ArrayView<byte> results)
|
||||
{
|
||||
ulong divisor = divisors[index];
|
||||
if (divisor <= 1UL)
|
||||
{
|
||||
results[index] = 0;
|
||||
return;
|
||||
}
|
||||
|
||||
results[index] = n % divisor == 0UL ? (byte)1 : (byte)0;
|
||||
}
|
||||
|
||||
public static void SharesFactorKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<byte> results)
|
||||
{
|
||||
|
||||
@@ -112,7 +112,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
resources.HitIndexBuffer,
|
||||
resources.Divisors,
|
||||
resources.Exponents,
|
||||
resources.FilteredDivisors,
|
||||
resources.DivisorData,
|
||||
resources.Offsets,
|
||||
resources.Counts,
|
||||
@@ -177,7 +176,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
MemoryBuffer1D<int, Stride1D.Dense> hitIndexBuffer,
|
||||
ulong[] divisors,
|
||||
ulong[] exponents,
|
||||
ulong[] filteredDivisors,
|
||||
GpuDivisorPartialData[] divisorData,
|
||||
int[] offsets,
|
||||
int[] counts,
|
||||
@@ -188,131 +186,36 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
{
|
||||
int batchCapacity = (int)divisorDataBuffer.Length;
|
||||
bool composite = false;
|
||||
bool processedAll = false;
|
||||
processedCount = 0UL;
|
||||
lastProcessed = 0UL;
|
||||
|
||||
int chunkCapacity = Math.Max(1, batchCapacity);
|
||||
var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
|
||||
|
||||
UInt128 twoP128 = (UInt128)prime << 1;
|
||||
UInt128 allowedMax128 = allowedMax;
|
||||
UInt128 firstDivisor128 = twoP128 + UInt128.One;
|
||||
bool invalidStride = twoP128 == UInt128.Zero;
|
||||
bool outOfRange = firstDivisor128 > allowedMax128;
|
||||
UInt128 numerator = allowedMax128 - UInt128.One;
|
||||
UInt128 maxK128 = invalidStride || outOfRange ? UInt128.Zero : numerator / twoP128;
|
||||
bool hasCandidates = !invalidStride && !outOfRange && maxK128 != UInt128.Zero;
|
||||
if (!hasCandidates)
|
||||
Span<ulong> divisorSpan = divisors.AsSpan();
|
||||
Span<ulong> exponentSpan = exponents.AsSpan();
|
||||
Span<GpuDivisorPartialData> divisorDataSpan = divisorData.AsSpan();
|
||||
Span<int> offsetSpan = offsets.AsSpan();
|
||||
Span<int> countSpan = counts.AsSpan();
|
||||
Span<ulong> cycleSpan = cycles.AsSpan();
|
||||
|
||||
while (!composite)
|
||||
{
|
||||
coveredRange = true;
|
||||
return false;
|
||||
}
|
||||
|
||||
ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
|
||||
ulong currentK = 1UL;
|
||||
|
||||
byte step10 = (byte)(twoP128 % 10UL);
|
||||
byte step8 = (byte)(twoP128 % 8UL);
|
||||
byte step5 = (byte)(twoP128 % 5UL);
|
||||
byte step3 = (byte)(twoP128 % 3UL);
|
||||
byte step7 = (byte)(twoP128 % 7UL);
|
||||
byte step11 = (byte)(twoP128 % 11UL);
|
||||
|
||||
UInt128 currentDivisor128 = firstDivisor128;
|
||||
byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL));
|
||||
byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL));
|
||||
byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL));
|
||||
byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL));
|
||||
byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL));
|
||||
byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL));
|
||||
bool lastIsSeven = (prime & 3UL) == 3UL;
|
||||
|
||||
Span<ulong> filteredStorage = filteredDivisors.AsSpan();
|
||||
Span<ulong> divisorStorage = divisors.AsSpan();
|
||||
Span<ulong> exponentStorage = exponents.AsSpan();
|
||||
Span<GpuDivisorPartialData> divisorDataStorage = divisorData.AsSpan();
|
||||
Span<int> offsetStorage = offsets.AsSpan();
|
||||
Span<int> countStorage = counts.AsSpan();
|
||||
Span<ulong> cycleStorage = cycles.AsSpan();
|
||||
while (currentK <= maxK && !composite)
|
||||
{
|
||||
int chunkCount = Math.Min(chunkCapacity, batchCapacity);
|
||||
ulong remainingK = maxK - currentK + 1UL;
|
||||
if ((ulong)chunkCount > remainingK)
|
||||
{
|
||||
chunkCount = (int)remainingK;
|
||||
}
|
||||
|
||||
if (chunkCount <= 0)
|
||||
{
|
||||
processedAll = true;
|
||||
break;
|
||||
}
|
||||
|
||||
int filteredCount = 0;
|
||||
UInt128 nextDivisor128 = currentDivisor128;
|
||||
|
||||
byte localRemainder10 = remainder10;
|
||||
byte localRemainder8 = remainder8;
|
||||
byte localRemainder5 = remainder5;
|
||||
byte localRemainder3 = remainder3;
|
||||
byte localRemainder7 = remainder7;
|
||||
byte localRemainder11 = remainder11;
|
||||
|
||||
for (int i = 0; i < chunkCount; i++)
|
||||
{
|
||||
ulong candidate = (ulong)nextDivisor128;
|
||||
|
||||
if (CandidatePassesHeuristics(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastIsSeven))
|
||||
{
|
||||
filteredStorage[filteredCount++] = candidate;
|
||||
}
|
||||
|
||||
nextDivisor128 += twoP128;
|
||||
localRemainder10 = AddMod(localRemainder10, step10, 10);
|
||||
localRemainder8 = AddMod(localRemainder8, step8, 8);
|
||||
localRemainder5 = AddMod(localRemainder5, step5, 5);
|
||||
localRemainder3 = AddMod(localRemainder3, step3, 3);
|
||||
localRemainder7 = AddMod(localRemainder7, step7, 7);
|
||||
localRemainder11 = AddMod(localRemainder11, step11, 11);
|
||||
}
|
||||
|
||||
remainder10 = localRemainder10;
|
||||
remainder8 = localRemainder8;
|
||||
remainder5 = localRemainder5;
|
||||
remainder3 = localRemainder3;
|
||||
remainder7 = localRemainder7;
|
||||
remainder11 = localRemainder11;
|
||||
|
||||
processedCount += (ulong)chunkCount;
|
||||
|
||||
UInt128 lastDivisor128 = nextDivisor128 - twoP128;
|
||||
lastProcessed = (ulong)lastDivisor128;
|
||||
currentDivisor128 = nextDivisor128;
|
||||
|
||||
currentK += (ulong)chunkCount;
|
||||
processedAll = currentK > maxK;
|
||||
|
||||
if (filteredCount == 0)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
Span<ulong> filteredDivisorsSpan = filteredStorage[..filteredCount];
|
||||
|
||||
Span<ulong> divisorSpan = divisorStorage;
|
||||
Span<ulong> exponentSpan = exponentStorage;
|
||||
Span<GpuDivisorPartialData> divisorDataSpan = divisorDataStorage;
|
||||
Span<int> offsetSpan = offsetStorage;
|
||||
Span<int> countSpan = countStorage;
|
||||
Span<ulong> cycleSpan = cycleStorage;
|
||||
int admissibleCount = 0;
|
||||
int exponentIndex = 0;
|
||||
for (int i = 0; i < filteredCount; i++)
|
||||
|
||||
while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
|
||||
{
|
||||
ulong divisorValue = filteredDivisorsSpan[i];
|
||||
MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
|
||||
ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData);
|
||||
processedCount = enumerator.ProcessedCount;
|
||||
ulong divisorValue = candidate.Value;
|
||||
|
||||
PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate);
|
||||
ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength(
|
||||
prime,
|
||||
in preparation,
|
||||
out _,
|
||||
out _,
|
||||
out _);
|
||||
|
||||
if (divisorCycle != prime)
|
||||
{
|
||||
continue;
|
||||
@@ -329,8 +232,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
exponentIndex++;
|
||||
}
|
||||
|
||||
processedCount = enumerator.ProcessedCount;
|
||||
lastProcessed = enumerator.LastDivisor;
|
||||
|
||||
if (admissibleCount == 0)
|
||||
{
|
||||
if (enumerator.Exhausted)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
continue;
|
||||
}
|
||||
|
||||
@@ -362,74 +273,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex];
|
||||
}
|
||||
|
||||
coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128);
|
||||
coveredRange = composite || enumerator.Exhausted;
|
||||
return composite;
|
||||
}
|
||||
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static ulong ResolveDivisorCycle(ulong divisor, ulong prime, in MontgomeryDivisorData divisorData)
|
||||
{
|
||||
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
|
||||
{
|
||||
return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
|
||||
}
|
||||
|
||||
return computedCycle;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static bool CandidatePassesHeuristics(
|
||||
byte remainder10,
|
||||
byte remainder8,
|
||||
byte remainder3,
|
||||
byte remainder5,
|
||||
byte remainder7,
|
||||
byte remainder11,
|
||||
bool lastIsSeven)
|
||||
{
|
||||
const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
|
||||
const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
|
||||
|
||||
ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
|
||||
if (((decimalMask >> remainder10) & 1) == 0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
if (remainder8 != 1 && remainder8 != 7)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static byte AddMod(byte value, byte delta, byte modulus)
|
||||
{
|
||||
int result = value + delta;
|
||||
if (result >= modulus)
|
||||
{
|
||||
result -= modulus;
|
||||
if (result >= modulus)
|
||||
{
|
||||
result -= modulus;
|
||||
}
|
||||
}
|
||||
|
||||
return (byte)result;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime)
|
||||
{
|
||||
// The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1,
|
||||
@@ -720,8 +567,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
|
||||
internal ulong[] Exponents { get; private set; } = null!;
|
||||
|
||||
internal ulong[] FilteredDivisors { get; private set; } = null!;
|
||||
|
||||
internal GpuDivisorPartialData[] DivisorData { get; private set; } = null!;
|
||||
|
||||
internal int[] Offsets { get; private set; } = null!;
|
||||
@@ -796,7 +641,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
|
||||
Divisors = ulongPool.Rent(capacity);
|
||||
Exponents = ulongPool.Rent(capacity);
|
||||
FilteredDivisors = ulongPool.Rent(capacity);
|
||||
DivisorData = partialPool.Rent(capacity);
|
||||
Offsets = intPool.Rent(capacity);
|
||||
Counts = intPool.Rent(capacity);
|
||||
@@ -835,14 +679,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
|
||||
ArrayPool<GpuDivisorPartialData> partialPool = ThreadStaticPools.GpuDivisorPartialDataPool;
|
||||
ulongPool.Return(Divisors, clearArray: false);
|
||||
ulongPool.Return(Exponents, clearArray: false);
|
||||
ulongPool.Return(FilteredDivisors, clearArray: false);
|
||||
partialPool.Return(DivisorData, clearArray: false);
|
||||
intPool.Return(Offsets, clearArray: false);
|
||||
intPool.Return(Counts, clearArray: false);
|
||||
ulongPool.Return(Cycles, clearArray: false);
|
||||
Divisors = Array.Empty<ulong>();
|
||||
Exponents = Array.Empty<ulong>();
|
||||
FilteredDivisors = Array.Empty<ulong>();
|
||||
DivisorData = Array.Empty<GpuDivisorPartialData>();
|
||||
Offsets = Array.Empty<int>();
|
||||
Counts = Array.Empty<int>();
|
||||
|
||||
@@ -0,0 +1,209 @@
|
||||
using System;
|
||||
|
||||
namespace PerfectNumbers.Core;
|
||||
|
||||
public static class HeuristicPrimeSieves
|
||||
{
|
||||
private const int SieveLength = PerfectNumberConstants.MaxQForDivisorCycles;
|
||||
private const int Wheel210 = 210;
|
||||
|
||||
private static readonly int[] GroupAConstants = [3, 7, 11, 13];
|
||||
private static readonly int[] GroupAIncrementPattern = [20, 10];
|
||||
|
||||
private static readonly ushort[] GroupBResiduesEnding1 = [1, 11, 31, 41, 61, 71, 101, 121, 131, 151, 181, 191];
|
||||
private static readonly ushort[] GroupBResiduesEnding3 = [13, 23, 43, 53, 73, 83, 103, 113, 143, 163, 173, 193];
|
||||
private static readonly ushort[] GroupBResiduesEnding7 = [17, 37, 47, 67, 97, 107, 127, 137, 157, 167, 187, 197];
|
||||
private static readonly ushort[] GroupBResiduesEnding9 = [19, 29, 59, 79, 89, 109, 139, 149, 169, 179, 199, 209];
|
||||
|
||||
private static int[]? s_groupADivisors;
|
||||
private static ulong[]? s_groupADivisorSquares;
|
||||
|
||||
private static int[]? s_groupBDivisorsEnding1;
|
||||
private static ulong[]? s_groupBDivisorSquaresEnding1;
|
||||
|
||||
private static int[]? s_groupBDivisorsEnding3;
|
||||
private static ulong[]? s_groupBDivisorSquaresEnding3;
|
||||
|
||||
private static int[]? s_groupBDivisorsEnding7;
|
||||
private static ulong[]? s_groupBDivisorSquaresEnding7;
|
||||
|
||||
private static int[]? s_groupBDivisorsEnding9;
|
||||
private static ulong[]? s_groupBDivisorSquaresEnding9;
|
||||
|
||||
private static readonly object InitLock = new();
|
||||
|
||||
public static ReadOnlySpan<int> GroupADivisors
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupAInitialized();
|
||||
return s_groupADivisors!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<ulong> GroupADivisorSquares
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupAInitialized();
|
||||
return s_groupADivisorSquares!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<int> GroupBDivisorsEnding1
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding1, ref s_groupBDivisorSquaresEnding1, GroupBResiduesEnding1);
|
||||
return s_groupBDivisorsEnding1!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding1
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding1, ref s_groupBDivisorSquaresEnding1, GroupBResiduesEnding1);
|
||||
return s_groupBDivisorSquaresEnding1!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<int> GroupBDivisorsEnding3
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding3, ref s_groupBDivisorSquaresEnding3, GroupBResiduesEnding3);
|
||||
return s_groupBDivisorsEnding3!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding3
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding3, ref s_groupBDivisorSquaresEnding3, GroupBResiduesEnding3);
|
||||
return s_groupBDivisorSquaresEnding3!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<int> GroupBDivisorsEnding7
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding7, ref s_groupBDivisorSquaresEnding7, GroupBResiduesEnding7);
|
||||
return s_groupBDivisorsEnding7!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding7
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding7, ref s_groupBDivisorSquaresEnding7, GroupBResiduesEnding7);
|
||||
return s_groupBDivisorSquaresEnding7!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<int> GroupBDivisorsEnding9
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding9, ref s_groupBDivisorSquaresEnding9, GroupBResiduesEnding9);
|
||||
return s_groupBDivisorsEnding9!;
|
||||
}
|
||||
}
|
||||
|
||||
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding9
|
||||
{
|
||||
get
|
||||
{
|
||||
EnsureGroupBInitialized(ref s_groupBDivisorsEnding9, ref s_groupBDivisorSquaresEnding9, GroupBResiduesEnding9);
|
||||
return s_groupBDivisorSquaresEnding9!;
|
||||
}
|
||||
}
|
||||
|
||||
private static void EnsureGroupAInitialized()
|
||||
{
|
||||
if (s_groupADivisors is not null)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
lock (InitLock)
|
||||
{
|
||||
if (s_groupADivisors is not null)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
var divisors = new int[SieveLength];
|
||||
var squares = new ulong[SieveLength];
|
||||
|
||||
int index = 0;
|
||||
for (; index < GroupAConstants.Length && index < divisors.Length; index++)
|
||||
{
|
||||
int value = GroupAConstants[index];
|
||||
divisors[index] = value;
|
||||
squares[index] = (ulong)value * (ulong)value;
|
||||
}
|
||||
|
||||
int candidate = 23;
|
||||
int incrementIndex = 0;
|
||||
while (index < divisors.Length)
|
||||
{
|
||||
divisors[index] = candidate;
|
||||
squares[index] = (ulong)candidate * (ulong)candidate;
|
||||
|
||||
candidate += GroupAIncrementPattern[incrementIndex];
|
||||
incrementIndex ^= 1;
|
||||
index++;
|
||||
}
|
||||
|
||||
s_groupADivisors = divisors;
|
||||
s_groupADivisorSquares = squares;
|
||||
}
|
||||
}
|
||||
|
||||
private static void EnsureGroupBInitialized(ref int[]? divisorsField, ref ulong[]? squaresField, ReadOnlySpan<ushort> residues)
|
||||
{
|
||||
if (divisorsField is not null)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
lock (InitLock)
|
||||
{
|
||||
if (divisorsField is not null)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
var divisors = new int[SieveLength];
|
||||
var squares = new ulong[SieveLength];
|
||||
|
||||
int index = 0;
|
||||
ulong baseValue = 0UL;
|
||||
|
||||
while (index < divisors.Length)
|
||||
{
|
||||
for (int i = 0; i < residues.Length && index < divisors.Length; i++)
|
||||
{
|
||||
ulong candidate = baseValue + residues[i];
|
||||
if (candidate <= 13UL)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
divisors[index] = (int)candidate;
|
||||
squares[index] = candidate * candidate;
|
||||
index++;
|
||||
}
|
||||
|
||||
baseValue += (ulong)Wheel210;
|
||||
}
|
||||
|
||||
divisorsField = divisors;
|
||||
squaresField = squares;
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,104 @@
|
||||
using System.Runtime.CompilerServices;
|
||||
using PerfectNumbers.Core.Gpu;
|
||||
|
||||
namespace PerfectNumbers.Core;
|
||||
|
||||
internal struct MersenneDivisorResidueStepper
|
||||
{
|
||||
private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
|
||||
private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
|
||||
|
||||
private readonly ushort decimalMask;
|
||||
private readonly byte step10;
|
||||
private readonly byte step8;
|
||||
private readonly byte step5;
|
||||
private readonly byte step3;
|
||||
private readonly byte step7;
|
||||
private readonly byte step11;
|
||||
|
||||
public byte Remainder10 { get; private set; }
|
||||
|
||||
public byte Remainder8 { get; private set; }
|
||||
|
||||
public byte Remainder5 { get; private set; }
|
||||
|
||||
public byte Remainder3 { get; private set; }
|
||||
|
||||
public byte Remainder7 { get; private set; }
|
||||
|
||||
public byte Remainder11 { get; private set; }
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public MersenneDivisorResidueStepper(ulong prime, in GpuUInt128 step, in GpuUInt128 firstDivisor)
|
||||
{
|
||||
decimalMask = (prime & 3UL) == 3UL ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
|
||||
|
||||
step10 = ComputeModulo(step, 10, multiplier: 6);
|
||||
step8 = (byte)(step.Low % 8UL);
|
||||
step5 = ComputeModulo(step, 5, multiplier: 1);
|
||||
step3 = ComputeModulo(step, 3, multiplier: 1);
|
||||
step7 = ComputeModulo(step, 7, multiplier: 2);
|
||||
step11 = ComputeModulo(step, 11, multiplier: 5);
|
||||
|
||||
Remainder10 = ComputeModulo(firstDivisor, 10, multiplier: 6);
|
||||
Remainder8 = (byte)(firstDivisor.Low % 8UL);
|
||||
Remainder5 = ComputeModulo(firstDivisor, 5, multiplier: 1);
|
||||
Remainder3 = ComputeModulo(firstDivisor, 3, multiplier: 1);
|
||||
Remainder7 = ComputeModulo(firstDivisor, 7, multiplier: 2);
|
||||
Remainder11 = ComputeModulo(firstDivisor, 11, multiplier: 5);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public readonly bool IsAdmissible()
|
||||
{
|
||||
if (((decimalMask >> Remainder10) & 1) == 0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
if (Remainder8 != 1 && Remainder8 != 7)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
if (Remainder3 == 0 || Remainder5 == 0 || Remainder7 == 0 || Remainder11 == 0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public void Advance()
|
||||
{
|
||||
Remainder10 = AddMod(Remainder10, step10, 10);
|
||||
Remainder8 = AddMod(Remainder8, step8, 8);
|
||||
Remainder5 = AddMod(Remainder5, step5, 5);
|
||||
Remainder3 = AddMod(Remainder3, step3, 3);
|
||||
Remainder7 = AddMod(Remainder7, step7, 7);
|
||||
Remainder11 = AddMod(Remainder11, step11, 11);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static byte ComputeModulo(in GpuUInt128 value, byte modulus, byte multiplier)
|
||||
{
|
||||
ulong high = value.High % modulus;
|
||||
ulong low = value.Low % modulus;
|
||||
int result = (int)((high * multiplier) + low);
|
||||
result %= modulus;
|
||||
return (byte)result;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static byte AddMod(byte value, byte delta, byte modulus)
|
||||
{
|
||||
int sum = value + delta;
|
||||
if (sum >= modulus)
|
||||
{
|
||||
sum -= modulus;
|
||||
}
|
||||
|
||||
return (byte)sum;
|
||||
}
|
||||
}
|
||||
@@ -1055,6 +1055,8 @@ internal static partial class PrimeOrderCalculator
|
||||
List<PendingEntry> pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2);
|
||||
Stack<ulong>? compositeStack = null;
|
||||
PartialFactorResult result;
|
||||
PrimeTester primeTester = PrimeTester.Exclusive;
|
||||
|
||||
if (!gpuFactored)
|
||||
{
|
||||
// We don't need to worry about leftovers, because we always use indexes within the calculated counts
|
||||
@@ -1120,7 +1122,7 @@ internal static partial class PrimeOrderCalculator
|
||||
// continue;
|
||||
// }
|
||||
|
||||
bool isPrime = PrimeTester.IsPrimeGpu(composite);
|
||||
bool isPrime = primeTester.HeuristicIsPrimeGpu(composite);
|
||||
// bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
|
||||
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
|
||||
|
||||
@@ -1182,9 +1184,10 @@ internal static partial class PrimeOrderCalculator
|
||||
|
||||
if (!entry.HasKnownPrimality)
|
||||
{
|
||||
bool isPrime = PrimeTester.IsPrimeGpu(composite);
|
||||
bool isPrime = primeTester.HeuristicIsPrimeGpu(composite);
|
||||
// bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
|
||||
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
|
||||
|
||||
entry = entry.WithPrimality(isPrime);
|
||||
pending[index] = entry;
|
||||
}
|
||||
@@ -1213,7 +1216,7 @@ internal static partial class PrimeOrderCalculator
|
||||
}
|
||||
else
|
||||
{
|
||||
cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
|
||||
cofactorIsPrime = primeTester.HeuristicIsPrimeGpu(cofactor);
|
||||
// cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None);
|
||||
// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
|
||||
}
|
||||
@@ -1743,7 +1746,9 @@ internal static partial class PrimeOrderCalculator
|
||||
return;
|
||||
}
|
||||
|
||||
bool isPrime = PrimeTester.IsPrimeGpu(value);
|
||||
PrimeTester primeTester = PrimeTester.Exclusive;
|
||||
|
||||
bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
|
||||
// bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None);
|
||||
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
|
||||
if (!knownComposite && isPrime)
|
||||
@@ -1755,7 +1760,7 @@ internal static partial class PrimeOrderCalculator
|
||||
ulong factor = PollardRhoStrict(value);
|
||||
ulong quotient = value / factor;
|
||||
|
||||
isPrime = PrimeTester.IsPrimeGpu(factor);
|
||||
isPrime = primeTester.HeuristicIsPrimeGpu(factor);
|
||||
// isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None);
|
||||
// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
|
||||
if (isPrime)
|
||||
@@ -1786,7 +1791,7 @@ internal static partial class PrimeOrderCalculator
|
||||
}
|
||||
else
|
||||
{
|
||||
isPrime = PrimeTester.IsPrimeGpu(quotient);
|
||||
isPrime = primeTester.HeuristicIsPrimeGpu(quotient);
|
||||
// isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None);
|
||||
// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
|
||||
}
|
||||
|
||||
@@ -1,5 +1,8 @@
|
||||
using System;
|
||||
using System.Buffers;
|
||||
using System.Collections.Concurrent;
|
||||
using Open.Numeric.Primes;
|
||||
using System.Collections.Generic;
|
||||
using System.Numerics;
|
||||
using System.Runtime.CompilerServices;
|
||||
using ILGPU;
|
||||
@@ -11,15 +14,96 @@ namespace PerfectNumbers.Core;
|
||||
|
||||
public sealed class PrimeTester(bool useInternal = false)
|
||||
{
|
||||
public bool IsPrime(ulong n, CancellationToken ct) =>
|
||||
// TODO: Eliminate this wrapper once callers can reach IsPrimeInternal directly; the extra
|
||||
// indirection shows up when we sieve millions of candidates per second.
|
||||
IsPrimeInternal(n, ct);
|
||||
// TODO: Don't use these tiny arrays. Hard-code the values / checks instead
|
||||
private static readonly ulong[] GroupAConstantDivisors = [3UL, 7UL, 11UL, 13UL];
|
||||
private static readonly byte[] GroupAIncrementPattern = [20, 10];
|
||||
private static readonly byte[] GroupBEndingOrderMod1 = [9, 1];
|
||||
private static readonly byte[] GroupBEndingOrderMod3 = [9, 7];
|
||||
private static readonly byte[] GroupBEndingOrderMod7 = [7, 1];
|
||||
private static readonly byte[] GroupBEndingOrderMod9 = [9, 7, 1];
|
||||
|
||||
public static bool IsPrimeGpu(ulong n) => new PrimeTester().IsPrimeGpu(n, CancellationToken.None);
|
||||
// TODO: Would it speed up things if we implement something like ExponentRemainderStepper,
|
||||
// MersenneDivisorResidueStepper, or CycleRemainderStepper work with these?
|
||||
private static readonly ushort[] Wheel210ResiduesEnding1 = [1, 11, 31, 41, 61, 71, 101, 121, 131, 151, 181, 191];
|
||||
private static readonly ushort[] Wheel210ResiduesEnding7 = [17, 37, 47, 67, 97, 107, 127, 137, 157, 167, 187, 197];
|
||||
private static readonly ushort[] Wheel210ResiduesEnding9 = [19, 29, 59, 79, 89, 109, 139, 149, 169, 179, 199, 209];
|
||||
private const ulong Wheel210 = 210UL;
|
||||
internal const int MaxGroupBSequences = 36;
|
||||
|
||||
private const bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B.
|
||||
|
||||
private readonly bool _useLegacyPrimeTester = useInternal;
|
||||
|
||||
private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
|
||||
|
||||
internal enum HeuristicDivisorGroup : byte
|
||||
{
|
||||
None = 0,
|
||||
GroupAConstant = 1,
|
||||
GroupAWheel = 2,
|
||||
GroupB = 3,
|
||||
}
|
||||
|
||||
internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue)
|
||||
{
|
||||
public readonly ulong Value = value;
|
||||
public readonly HeuristicDivisorGroup Group = group;
|
||||
public readonly byte Ending = ending;
|
||||
public readonly byte PriorityIndex = priorityIndex;
|
||||
public readonly ushort WheelResidue = wheelResidue;
|
||||
}
|
||||
|
||||
internal readonly struct HeuristicDivisorPreparation(
|
||||
in HeuristicDivisorCandidate candidate,
|
||||
in MontgomeryDivisorData divisorData,
|
||||
ulong cycleLengthHint,
|
||||
bool hasCycleLengthHint)
|
||||
{
|
||||
public readonly HeuristicDivisorCandidate Candidate = candidate;
|
||||
public readonly MontgomeryDivisorData DivisorData = divisorData;
|
||||
public readonly ulong CycleLengthHint = cycleLengthHint;
|
||||
public readonly bool HasCycleLengthHint = hasCycleLengthHint;
|
||||
public readonly bool RequiresCycleComputation = !hasCycleLengthHint;
|
||||
}
|
||||
|
||||
public bool IsPrime(ulong n, CancellationToken ct)
|
||||
{
|
||||
if (_useLegacyPrimeTester)
|
||||
{
|
||||
return LegacyIsPrimeInternal(n, ct);
|
||||
}
|
||||
|
||||
return IsPrimeInternal(n, ct);
|
||||
}
|
||||
|
||||
[ThreadStatic]
|
||||
private static PrimeTester? _tester;
|
||||
|
||||
public static PrimeTester Exclusive => _tester ??= new();
|
||||
|
||||
public static bool IsPrimeGpu(ulong n)
|
||||
{
|
||||
return Exclusive.IsPrimeGpu(n, CancellationToken.None);
|
||||
}
|
||||
|
||||
public static bool IsPrimeGpu(ulong n, ulong limit, byte nMod10)
|
||||
{
|
||||
return Exclusive.HeuristicIsPrimeGpu(n, limit, nMod10);
|
||||
}
|
||||
|
||||
// Optional GPU-assisted primality: batched small-prime sieve on device.
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public bool IsPrimeGpu(ulong n, CancellationToken ct)
|
||||
{
|
||||
if (_useLegacyPrimeTester || !EnableHeuristicPrimeTesting)
|
||||
{
|
||||
return LegacyIsPrimeGpu(n, ct);
|
||||
}
|
||||
|
||||
return HeuristicIsPrimeGpu(n);
|
||||
}
|
||||
|
||||
public static bool LegacyIsPrimeGpu(ulong n, CancellationToken ct)
|
||||
{
|
||||
// EvenPerfectBitScanner never enqueues exponents <= 3 on the GPU path; keep the legacy guard commented out
|
||||
// so ad-hoc callers that expect that behavior must re-enable it intentionally.
|
||||
@@ -34,11 +118,11 @@ public sealed class PrimeTester(bool useInternal = false)
|
||||
|
||||
// EvenPerfectBitScanner seeds p at 136,279,841 and advances monotonically, so production scans never hit the GPU with
|
||||
// exponents below 31. Leave this CPU redirect commented out to keep the wrapper branchless while documenting the guard.
|
||||
// The ternary fallback below still routes sub-31 values through IsPrimeInternal for correctness when ad-hoc callers bypass
|
||||
// The ternary fallback below still routes sub-31 values through LegacyIsPrimeInternal for correctness when ad-hoc callers bypass
|
||||
// the scanner pipeline.
|
||||
// if (n < 31UL)
|
||||
// {
|
||||
// return IsPrimeInternal(n, ct);
|
||||
// return LegacyIsPrimeInternal(n, ct);
|
||||
// }
|
||||
|
||||
bool forceCpu = GpuContextPool.ForceCpu;
|
||||
@@ -58,16 +142,611 @@ public sealed class PrimeTester(bool useInternal = false)
|
||||
bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0;
|
||||
bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime;
|
||||
|
||||
// Defensive fallback: GPU sieve produced a composite verdict or execution was skipped.
|
||||
// Confirm with CPU logic to prevent false negatives observed on some accelerators.
|
||||
return requiresCpuFallback ? IsPrimeInternal(n, ct) : true;
|
||||
return requiresCpuFallback ? LegacyIsPrimeInternal(n, ct) : true;
|
||||
}
|
||||
|
||||
public bool HeuristicIsPrimeGpu(ulong n)
|
||||
{
|
||||
byte nMod10 = (byte)n.Mod10();
|
||||
ulong sqrtLimit = ComputeHeuristicSqrt(n);
|
||||
return HeuristicIsPrimeGpu(n, sqrtLimit, nMod10);
|
||||
}
|
||||
|
||||
private bool HeuristicIsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10)
|
||||
{
|
||||
// TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
|
||||
if (sqrtLimit < 3UL)
|
||||
{
|
||||
return EvaluateWithOpenNumericFallback(n);
|
||||
}
|
||||
|
||||
if (!UseHeuristicGroupBTrialDivision)
|
||||
{
|
||||
return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false);
|
||||
}
|
||||
|
||||
if (GpuContextPool.ForceCpu)
|
||||
{
|
||||
return HeuristicIsPrimeCpu(n, sqrtLimit, nMod10);
|
||||
}
|
||||
|
||||
bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10);
|
||||
return !compositeDetected;
|
||||
}
|
||||
|
||||
internal static bool HeuristicIsPrimeCpu(ulong n)
|
||||
{
|
||||
return HeuristicIsPrimeCpu(n, 0UL, (byte)n.Mod10());
|
||||
}
|
||||
|
||||
internal static bool HeuristicIsPrimeCpu(ulong n, ulong sqrtLimit, byte nMod10)
|
||||
{
|
||||
if (sqrtLimit == 0UL)
|
||||
{
|
||||
sqrtLimit = ComputeHeuristicSqrt(n);
|
||||
}
|
||||
|
||||
if (sqrtLimit < 3UL)
|
||||
{
|
||||
return EvaluateWithOpenNumericFallback(n);
|
||||
}
|
||||
|
||||
bool includeGroupB = UseHeuristicGroupBTrialDivision;
|
||||
return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB);
|
||||
}
|
||||
|
||||
private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB)
|
||||
{
|
||||
Span<HeuristicGroupBSequenceState> groupBBuffer = stackalloc HeuristicGroupBSequenceState[MaxGroupBSequences];
|
||||
var enumerator = new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
|
||||
|
||||
while (enumerator.TryGetNext(out HeuristicDivisorCandidate candidate))
|
||||
{
|
||||
if (!includeGroupB && candidate.Group == HeuristicDivisorGroup.GroupB)
|
||||
{
|
||||
return EvaluateWithOpenNumericFallback(n);
|
||||
}
|
||||
|
||||
ulong divisor = candidate.Value;
|
||||
if (divisor <= 1UL)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
if (n % divisor == 0UL)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
return includeGroupB ? true : EvaluateWithOpenNumericFallback(n);
|
||||
}
|
||||
|
||||
private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10)
|
||||
{
|
||||
Span<HeuristicGroupBSequenceState> groupBBuffer = stackalloc HeuristicGroupBSequenceState[MaxGroupBSequences];
|
||||
var enumerator = new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
|
||||
|
||||
int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize);
|
||||
var candidatePool = ArrayPool<HeuristicDivisorCandidate>.Shared;
|
||||
var divisorPool = ThreadStaticPools.UlongPool;
|
||||
var hitPool = ThreadStaticPools.BytePool;
|
||||
|
||||
HeuristicDivisorCandidate[]? candidateArray = null;
|
||||
ulong[]? divisorArray = null;
|
||||
byte[]? hitFlags = null;
|
||||
|
||||
var limiter = GpuPrimeWorkLimiter.Acquire();
|
||||
var gpu = PrimeTesterGpuContextPool.Rent();
|
||||
|
||||
bool compositeDetected = false;
|
||||
|
||||
try
|
||||
{
|
||||
candidateArray = candidatePool.Rent(batchCapacity);
|
||||
divisorArray = divisorPool.Rent(batchCapacity);
|
||||
hitFlags = hitPool.Rent(batchCapacity);
|
||||
|
||||
var accelerator = gpu.Accelerator;
|
||||
var state = GpuKernelState.GetOrCreate(accelerator);
|
||||
|
||||
lock (gpu.ExecutionLock)
|
||||
{
|
||||
var scratch = state.RentScratch(batchCapacity, accelerator);
|
||||
try
|
||||
{
|
||||
int count = 0;
|
||||
|
||||
bool ProcessBatch(int length)
|
||||
{
|
||||
scratch.Input.View.CopyFromCPU(ref divisorArray![0], length);
|
||||
state.HeuristicTrialDivisionKernel(length, scratch.Input.View, n, scratch.Output.View);
|
||||
accelerator.Synchronize();
|
||||
scratch.Output.View.CopyToCPU(ref hitFlags![0], length);
|
||||
|
||||
for (int i = 0; i < length; i++)
|
||||
{
|
||||
if (hitFlags![i] == 0)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
ulong divisor = candidateArray![i].Value;
|
||||
if (divisor > 1UL && n % divisor == 0UL)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
while (enumerator.TryGetNext(out HeuristicDivisorCandidate candidate))
|
||||
{
|
||||
ulong divisor = candidate.Value;
|
||||
if (divisor <= 1UL)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
candidateArray[count] = candidate;
|
||||
divisorArray[count] = divisor;
|
||||
count++;
|
||||
|
||||
if (count == batchCapacity)
|
||||
{
|
||||
if (ProcessBatch(count))
|
||||
{
|
||||
compositeDetected = true;
|
||||
break;
|
||||
}
|
||||
|
||||
count = 0;
|
||||
}
|
||||
}
|
||||
|
||||
if (!compositeDetected && count > 0)
|
||||
{
|
||||
if (ProcessBatch(count))
|
||||
{
|
||||
compositeDetected = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
finally
|
||||
{
|
||||
state.ReturnScratch(scratch);
|
||||
}
|
||||
}
|
||||
}
|
||||
finally
|
||||
{
|
||||
if (hitFlags is not null)
|
||||
{
|
||||
hitPool.Return(hitFlags);
|
||||
}
|
||||
|
||||
if (divisorArray is not null)
|
||||
{
|
||||
divisorPool.Return(divisorArray);
|
||||
}
|
||||
|
||||
if (candidateArray is not null)
|
||||
{
|
||||
candidatePool.Return(candidateArray);
|
||||
}
|
||||
|
||||
gpu.Dispose();
|
||||
limiter.Dispose();
|
||||
}
|
||||
|
||||
return compositeDetected;
|
||||
}
|
||||
|
||||
private static bool EvaluateWithOpenNumericFallback(ulong n)
|
||||
{
|
||||
return Prime.Numbers.IsPrime(n);
|
||||
}
|
||||
|
||||
internal static ulong ComputeHeuristicSqrt(ulong n)
|
||||
{
|
||||
ulong sqrt = (ulong)Math.Sqrt(n);
|
||||
UInt128 square = (UInt128)sqrt * sqrt;
|
||||
|
||||
while (square > n)
|
||||
{
|
||||
sqrt--;
|
||||
square = (UInt128)sqrt * sqrt;
|
||||
}
|
||||
|
||||
ulong next = sqrt + 1UL;
|
||||
UInt128 nextSquare = (UInt128)next * next;
|
||||
while (nextSquare <= n)
|
||||
{
|
||||
sqrt = next;
|
||||
next++;
|
||||
nextSquare = (UInt128)next * next;
|
||||
}
|
||||
|
||||
return sqrt;
|
||||
}
|
||||
|
||||
private static ReadOnlySpan<byte> GetGroupBEndingOrder(byte nMod10) => nMod10 switch
|
||||
{
|
||||
1 => GroupBEndingOrderMod1,
|
||||
3 => GroupBEndingOrderMod3,
|
||||
7 => GroupBEndingOrderMod7,
|
||||
9 => GroupBEndingOrderMod9,
|
||||
_ => ReadOnlySpan<byte>.Empty,
|
||||
};
|
||||
|
||||
private static ReadOnlySpan<ushort> GetWheel210ResiduesForEnding(byte ending) => ending switch
|
||||
{
|
||||
1 => Wheel210ResiduesEnding1,
|
||||
7 => Wheel210ResiduesEnding7,
|
||||
9 => Wheel210ResiduesEnding9,
|
||||
_ => ReadOnlySpan<ushort>.Empty,
|
||||
};
|
||||
|
||||
private static int InitializeGroupBStates(byte nMod10, Span<HeuristicGroupBSequenceState> buffer)
|
||||
{
|
||||
ReadOnlySpan<byte> endings = GetGroupBEndingOrder(nMod10);
|
||||
int count = 0;
|
||||
|
||||
for (int priorityIndex = 0; priorityIndex < endings.Length; priorityIndex++)
|
||||
{
|
||||
ReadOnlySpan<ushort> residues = GetWheel210ResiduesForEnding(endings[priorityIndex]);
|
||||
for (int residueIndex = 0; residueIndex < residues.Length; residueIndex++)
|
||||
{
|
||||
buffer[count++] = new HeuristicGroupBSequenceState(residues[residueIndex], endings[priorityIndex], (byte)priorityIndex);
|
||||
}
|
||||
}
|
||||
|
||||
return count;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
|
||||
{
|
||||
return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor)
|
||||
{
|
||||
return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate)
|
||||
{
|
||||
ulong divisor = candidate.Value;
|
||||
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
|
||||
bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength);
|
||||
return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
internal static ulong ResolveHeuristicCycleLength(
|
||||
ulong exponent,
|
||||
in HeuristicDivisorPreparation preparation,
|
||||
out bool cycleFromHint,
|
||||
out bool cycleComputed,
|
||||
out bool primeOrderFailed)
|
||||
{
|
||||
if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL)
|
||||
{
|
||||
cycleFromHint = true;
|
||||
cycleComputed = true;
|
||||
primeOrderFailed = false;
|
||||
return preparation.CycleLengthHint;
|
||||
}
|
||||
|
||||
ulong divisor = preparation.Candidate.Value;
|
||||
MontgomeryDivisorData divisorData = preparation.DivisorData;
|
||||
|
||||
bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
|
||||
divisor,
|
||||
exponent,
|
||||
divisorData,
|
||||
out ulong computedCycle,
|
||||
out bool primeOrderFailedLocal);
|
||||
|
||||
if (trySuccess && computedCycle != 0UL)
|
||||
{
|
||||
cycleFromHint = false;
|
||||
cycleComputed = true;
|
||||
primeOrderFailed = primeOrderFailedLocal;
|
||||
return computedCycle;
|
||||
}
|
||||
|
||||
primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL;
|
||||
ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
|
||||
cycleFromHint = false;
|
||||
cycleComputed = resolvedCycle != 0UL;
|
||||
return resolvedCycle;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength)
|
||||
{
|
||||
ulong[] snapshot = HeuristicSmallCycleSnapshot;
|
||||
if (divisor < (ulong)snapshot.Length)
|
||||
{
|
||||
cycleLength = snapshot[(int)divisor];
|
||||
if (cycleLength != 0UL)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
cycleLength = 0UL;
|
||||
return false;
|
||||
}
|
||||
|
||||
internal ref struct HeuristicDivisorEnumerator
|
||||
{
|
||||
private readonly ulong sqrtLimit;
|
||||
private int groupAConstantIndex;
|
||||
private ulong groupAWheelCandidate;
|
||||
private int groupAIncrementIndex;
|
||||
private bool groupAWheelActive;
|
||||
private Span<HeuristicGroupBSequenceState> groupBStates;
|
||||
|
||||
public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
|
||||
{
|
||||
this.sqrtLimit = sqrtLimit;
|
||||
groupAConstantIndex = 0;
|
||||
groupAWheelCandidate = 23UL;
|
||||
groupAIncrementIndex = 0;
|
||||
groupAWheelActive = sqrtLimit >= 23UL;
|
||||
|
||||
int count = InitializeGroupBStates(nMod10, groupBBuffer);
|
||||
groupBStates = count == 0 ? Span<HeuristicGroupBSequenceState>.Empty : groupBBuffer[..count];
|
||||
|
||||
for (int i = 0; i < groupBStates.Length; i++)
|
||||
{
|
||||
groupBStates[i].TryNormalize();
|
||||
}
|
||||
}
|
||||
|
||||
public bool TryGetNext(out HeuristicDivisorCandidate candidate)
|
||||
{
|
||||
// TODO: Hard-code constant A divisor checks here instead of using an array
|
||||
while (groupAConstantIndex < GroupAConstantDivisors.Length)
|
||||
{
|
||||
ulong value = GroupAConstantDivisors[groupAConstantIndex++];
|
||||
if (value > sqrtLimit)
|
||||
{
|
||||
groupAConstantIndex = GroupAConstantDivisors.Length;
|
||||
break;
|
||||
}
|
||||
|
||||
// TODO: We're constantly dividing constants mod 10 which is constant too.
|
||||
candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAConstant, (byte)(value % 10UL), 0, (ushort)value);
|
||||
return true;
|
||||
}
|
||||
|
||||
if (groupAWheelActive)
|
||||
{
|
||||
ulong value = groupAWheelCandidate;
|
||||
if (value <= sqrtLimit)
|
||||
{
|
||||
candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAWheel, 3, 0, (ushort)(value % Wheel210));
|
||||
byte increment = GroupAIncrementPattern[groupAIncrementIndex];
|
||||
groupAIncrementIndex ^= 1;
|
||||
|
||||
ulong next = value + increment;
|
||||
if (next <= value)
|
||||
{
|
||||
groupAWheelActive = false;
|
||||
}
|
||||
else
|
||||
{
|
||||
groupAWheelCandidate = next;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
groupAWheelActive = false;
|
||||
}
|
||||
|
||||
ulong bestCandidate = ulong.MaxValue;
|
||||
int bestIndex = -1;
|
||||
|
||||
for (int i = 0; i < groupBStates.Length; i++)
|
||||
{
|
||||
ref HeuristicGroupBSequenceState state = ref groupBStates[i];
|
||||
ulong value = state.Candidate;
|
||||
if (value == ulong.MaxValue)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
if (value > sqrtLimit)
|
||||
{
|
||||
state.Candidate = ulong.MaxValue;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (bestIndex == -1 || value < bestCandidate ||
|
||||
(value == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex))
|
||||
{
|
||||
bestCandidate = value;
|
||||
bestIndex = i;
|
||||
}
|
||||
}
|
||||
|
||||
if (bestIndex == -1)
|
||||
{
|
||||
candidate = default;
|
||||
return false;
|
||||
}
|
||||
|
||||
ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex];
|
||||
var result = new HeuristicDivisorCandidate(bestState.Candidate, HeuristicDivisorGroup.GroupB, bestState.Ending, bestState.PriorityIndex, bestState.Residue);
|
||||
bestState.Advance();
|
||||
candidate = result;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
internal struct HeuristicGroupBSequenceState(ushort residue, byte ending, byte priorityIndex)
|
||||
{
|
||||
public ulong Candidate = residue;
|
||||
public byte PriorityIndex = priorityIndex;
|
||||
public byte Ending = ending;
|
||||
public ushort Residue = residue;
|
||||
|
||||
public void TryNormalize()
|
||||
{
|
||||
ulong candidate = Candidate;
|
||||
if (candidate == ulong.MaxValue)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
while (candidate <= 13UL)
|
||||
{
|
||||
ulong next = candidate + Wheel210;
|
||||
if (next <= candidate)
|
||||
{
|
||||
Candidate = ulong.MaxValue;
|
||||
return;
|
||||
}
|
||||
|
||||
candidate = next;
|
||||
}
|
||||
|
||||
Candidate = candidate;
|
||||
}
|
||||
|
||||
public void Advance()
|
||||
{
|
||||
ulong candidate = Candidate;
|
||||
if (candidate == ulong.MaxValue)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
ulong next = candidate + Wheel210;
|
||||
Candidate = next > candidate ? next : ulong.MaxValue;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
internal ref struct MersenneHeuristicDivisorEnumerator
|
||||
{
|
||||
private readonly GpuUInt128 step;
|
||||
private readonly GpuUInt128 limit;
|
||||
private GpuUInt128 current;
|
||||
private MersenneDivisorResidueStepper residueStepper;
|
||||
private bool active;
|
||||
private ulong processedCount;
|
||||
private ulong lastDivisor;
|
||||
|
||||
public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor)
|
||||
{
|
||||
var stepLocal = new GpuUInt128(exponent);
|
||||
stepLocal.ShiftLeft(1);
|
||||
step = stepLocal;
|
||||
|
||||
var limitLocal = new GpuUInt128(maxDivisor);
|
||||
limit = limitLocal;
|
||||
|
||||
var firstDivisor = stepLocal;
|
||||
firstDivisor.Add(1UL);
|
||||
|
||||
bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0;
|
||||
|
||||
current = hasCandidates ? firstDivisor : GpuUInt128.Zero;
|
||||
residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default;
|
||||
active = hasCandidates;
|
||||
processedCount = 0UL;
|
||||
lastDivisor = 0UL;
|
||||
}
|
||||
|
||||
public bool TryGetNext(out HeuristicDivisorCandidate candidate)
|
||||
{
|
||||
while (active)
|
||||
{
|
||||
ulong value = current.Low;
|
||||
processedCount++;
|
||||
lastDivisor = value;
|
||||
|
||||
bool admissible = residueStepper.IsAdmissible();
|
||||
|
||||
Advance();
|
||||
|
||||
if (admissible)
|
||||
{
|
||||
candidate = CreateCandidate(value);
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
candidate = default;
|
||||
return false;
|
||||
}
|
||||
|
||||
private void Advance()
|
||||
{
|
||||
GpuUInt128 next = current;
|
||||
next.Add(step);
|
||||
|
||||
if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0)
|
||||
{
|
||||
active = false;
|
||||
return;
|
||||
}
|
||||
|
||||
current = next;
|
||||
residueStepper.Advance();
|
||||
}
|
||||
|
||||
public ulong ProcessedCount => processedCount;
|
||||
|
||||
public ulong LastDivisor => lastDivisor;
|
||||
|
||||
public bool Exhausted => !active;
|
||||
|
||||
private static HeuristicDivisorCandidate CreateCandidate(ulong value)
|
||||
{
|
||||
byte ending = (byte)(value % 10UL);
|
||||
HeuristicDivisorGroup group;
|
||||
|
||||
if (value == 3UL || value == 7UL || value == 11UL || value == 13UL)
|
||||
{
|
||||
group = HeuristicDivisorGroup.GroupAConstant;
|
||||
}
|
||||
else if (ending == 3)
|
||||
{
|
||||
group = HeuristicDivisorGroup.GroupAWheel;
|
||||
}
|
||||
else
|
||||
{
|
||||
group = HeuristicDivisorGroup.GroupB;
|
||||
}
|
||||
|
||||
ushort residue = (ushort)(value % Wheel210);
|
||||
return new HeuristicDivisorCandidate(value, group, ending, 0, residue);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Note: GPU-backed primality path is implemented via IsPrimeGpu/IsPrimeBatchGpu and is routed
|
||||
// from EvenPerfectBitScanner based on --primes-device.
|
||||
|
||||
internal static bool IsPrimeInternal(ulong n, CancellationToken ct)
|
||||
{
|
||||
if (!EnableHeuristicPrimeTesting)
|
||||
{
|
||||
return LegacyIsPrimeInternal(n, ct);
|
||||
}
|
||||
|
||||
return HeuristicIsPrimeCpu(n);
|
||||
}
|
||||
|
||||
public static bool LegacyIsPrimeInternal(ulong n, CancellationToken ct)
|
||||
{
|
||||
// EvenPerfectBitScanner streams monotonically increasing odd exponents that start at 136,279,841 and already exclude
|
||||
// multiples of five. Factorization helpers (PrimeOrderCalculator, Pollard routines, and divisor-cycle warmups) reuse
|
||||
@@ -113,8 +792,13 @@ public sealed class PrimeTester(bool useInternal = false)
|
||||
return result;
|
||||
}
|
||||
|
||||
public static bool EnableHeuristicPrimeTesting { get; set; } = true;
|
||||
|
||||
public static int GpuBatchSize { get; set; } = 262_144;
|
||||
|
||||
public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096;
|
||||
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static void IsPrimeBatchGpu(ReadOnlySpan<ulong> values, Span<byte> results)
|
||||
{
|
||||
@@ -150,7 +834,7 @@ public sealed class PrimeTester(bool useInternal = false)
|
||||
int count = remaining > batchSize ? batchSize : remaining;
|
||||
|
||||
values.Slice(pos, count).CopyTo(temp);
|
||||
input.View.CopyFromCPU(ref temp[0], count);
|
||||
input.View.CopyFromCPU(ref temp[0], count);
|
||||
|
||||
state.Kernel(count, input.View, state.DevicePrimes.View, output.View);
|
||||
accelerator.Synchronize();
|
||||
@@ -263,6 +947,7 @@ public sealed class PrimeTester(bool useInternal = false)
|
||||
private sealed class KernelState
|
||||
{
|
||||
public Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>> Kernel { get; }
|
||||
public Action<Index1D, ArrayView<ulong>, ulong, ArrayView<byte>> HeuristicTrialDivisionKernel { get; }
|
||||
public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimes { get; }
|
||||
private readonly Accelerator _accel;
|
||||
private readonly System.Collections.Concurrent.ConcurrentBag<ScratchBuffers> _scratchPool = [];
|
||||
@@ -275,6 +960,7 @@ public sealed class PrimeTester(bool useInternal = false)
|
||||
_accel = accelerator;
|
||||
// Compile once per accelerator and upload primes once.
|
||||
Kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
|
||||
HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ulong, ArrayView<byte>>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
|
||||
|
||||
var primes = PrimesGenerator.SmallPrimes;
|
||||
DevicePrimes = accelerator.Allocate1D<uint>(primes.Length);
|
||||
|
||||
Reference in new issue
Block a user