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

Align heuristic prime testers with squared limits

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-10-29 19:35:59 +01:00
1 parent 5536dc4865
commit f7a0627ce5
4 files changed
+109 -139
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs
+4 -8
View File
@@ -19,7 +19,7 @@ public class HeuristicCombinedPrimeTesterTests
Span<HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState> groupBState = Span<HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState> groupBState =
stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4]; stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4];
var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500, lastDigit, groupBState); var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500UL * 500UL, lastDigit, groupBState);
var results = new List<HeuristicCombinedPrimeTester.HeuristicDivisorCandidate>(32); var results = new List<HeuristicCombinedPrimeTester.HeuristicDivisorCandidate>(32);
while (results.Count < 20 && enumerator.TryGetNext(out var candidate)) while (results.Count < 20 && enumerator.TryGetNext(out var candidate))
@@ -230,14 +230,10 @@ public class HeuristicCombinedPrimeTesterTests
[InlineData(18_446_744_073_709_551_615UL)] [InlineData(18_446_744_073_709_551_615UL)]
[InlineData(18_446_744_073_709_551_614UL)] [InlineData(18_446_744_073_709_551_614UL)]
[InlineData(18_446_744_073_709_551_613UL)] [InlineData(18_446_744_073_709_551_613UL)]
public void ComputeHeuristicSqrt_returns_floor_square_root(ulong value) public void ComputeHeuristicDivisorSquareLimit_returns_input_value(ulong value)
{ {
ulong sqrt = HeuristicCombinedPrimeTester.ComputeHeuristicSqrt(value); ulong limit = HeuristicCombinedPrimeTester.ComputeHeuristicDivisorSquareLimit(value);
bool squareBelowOrEqual = (UInt128)sqrt * sqrt <= (UInt128)value; limit.Should().Be(value);
bool nextSquareAbove = (UInt128)(sqrt + 1UL) * (UInt128)(sqrt + 1UL) > (UInt128)value;
squareBelowOrEqual.Should().BeTrue();
nextSquareAbove.Should().BeTrue();
} }
} }
PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs
+57 -70
View File
@@ -152,36 +152,36 @@ public sealed class HeuristicCombinedPrimeTester
public static bool IsPrimeCpu(ulong n) public static bool IsPrimeCpu(ulong n)
{ {
byte nMod10 = (byte)n.Mod10(); byte nMod10 = (byte)n.Mod10();
ulong sqrtLimit = ComputeHeuristicSqrt(n); ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
if (sqrtLimit < 3UL) if (maxDivisorSquare < 9UL)
{ {
return EvaluateWithOpenNumericFallback(n); return EvaluateWithOpenNumericFallback(n);
} }
return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10);
} }
public static bool IsPrimeGpu(ulong n) public static bool IsPrimeGpu(ulong n)
{ {
byte nMod10 = (byte)n.Mod10(); byte nMod10 = (byte)n.Mod10();
ulong sqrtLimit = ComputeHeuristicSqrt(n); ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
return IsPrimeGpu(n, sqrtLimit, nMod10); return IsPrimeGpu(n, maxDivisorSquare, nMod10);
} }
public static bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) public static bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10)
{ {
// TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
if (sqrtLimit < 3UL) if (maxDivisorSquare < 9UL)
{ {
return EvaluateWithOpenNumericFallback(n); return EvaluateWithOpenNumericFallback(n);
} }
bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10);
return !compositeDetected; return !compositeDetected;
} }
private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10)
{ {
ReadOnlySpan<uint> combinedDivisors = GetCombinedDivisors(nMod10); ReadOnlySpan<uint> combinedDivisors = GetCombinedDivisors(nMod10);
ReadOnlySpan<ulong> combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); ReadOnlySpan<ulong> combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
@@ -191,7 +191,7 @@ public sealed class HeuristicCombinedPrimeTester
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
ulong divisorSquare = combinedDivisorSquares[i]; ulong divisorSquare = combinedDivisorSquares[i];
if (divisorSquare > n) if (divisorSquare > maxDivisorSquare)
{ {
break; break;
} }
@@ -205,7 +205,7 @@ public sealed class HeuristicCombinedPrimeTester
return EvaluateWithOpenNumericFallback(n); return EvaluateWithOpenNumericFallback(n);
} }
private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
{ {
int batchCapacity = HeuristicGpuDivisorBatchSize; int batchCapacity = HeuristicGpuDivisorBatchSize;
var divisorPool = ThreadStaticPools.UlongPool; var divisorPool = ThreadStaticPools.UlongPool;
@@ -243,7 +243,7 @@ public sealed class HeuristicCombinedPrimeTester
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
ulong divisorSquare = combinedDivisorSquares[i]; ulong divisorSquare = combinedDivisorSquares[i];
if (divisorSquare > n) if (divisorSquare > maxDivisorSquare)
{ {
break; break;
} }
@@ -282,29 +282,10 @@ Cleanup:
{ {
return Prime.Numbers.IsPrime(n); return Prime.Numbers.IsPrime(n);
} }
internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n)
internal static ulong ComputeHeuristicSqrt(ulong n) {
{ return 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<uint> GetCombinedDivisors(byte nMod10) private static ReadOnlySpan<uint> GetCombinedDivisors(byte nMod10)
{ {
@@ -514,9 +495,9 @@ Cleanup:
}; };
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer) internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
{ {
return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer);
} }
internal struct HeuristicGroupBSequenceState internal struct HeuristicGroupBSequenceState
@@ -597,44 +578,50 @@ Cleanup:
} }
internal ref struct HeuristicDivisorEnumerator internal ref struct HeuristicDivisorEnumerator
{ {
private readonly ulong sqrtLimit; private readonly ulong maxDivisorSquare;
private readonly ReadOnlySpan<uint> combinedDivisors; private readonly ReadOnlySpan<uint> combinedDivisors;
private int index; private readonly ReadOnlySpan<ulong> combinedDivisorSquares;
private int index;
public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer) public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
{ {
this.sqrtLimit = sqrtLimit; this.maxDivisorSquare = maxDivisorSquare;
_ = groupBBuffer; _ = groupBBuffer;
combinedDivisors = GetCombinedDivisors(nMod10); combinedDivisors = GetCombinedDivisors(nMod10);
index = 0; combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
} index = 0;
}
public bool TryGetNext(out HeuristicDivisorCandidate candidate) public bool TryGetNext(out HeuristicDivisorCandidate candidate)
{ {
while (index < combinedDivisors.Length) while (index < combinedDivisors.Length)
{ {
uint entry = combinedDivisors[index++]; int currentIndex = index;
if (entry > sqrtLimit) ulong divisorSquare = combinedDivisorSquares[currentIndex];
{ if (divisorSquare > maxDivisorSquare)
index = combinedDivisors.Length; {
break; index = combinedDivisors.Length;
} break;
}
HeuristicDivisorGroup group = ResolveGroup(entry); uint entry = combinedDivisors[currentIndex];
ulong value = entry; index = currentIndex + 1;
byte ending = (byte)(value % 10UL);
ushort residue = (ushort)(value % Wheel210);
candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue);
return true;
}
candidate = default; HeuristicDivisorGroup group = ResolveGroup(entry);
return false; ulong value = entry;
} byte ending = (byte)(value % 10UL);
} ushort residue = (ushort)(value % Wheel210);
candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue);
return true;
}
internal ref struct MersenneHeuristicDivisorEnumerator candidate = default;
return false;
}
}
internal ref struct MersenneHeuristicDivisorEnumerator
{ {
private readonly GpuUInt128 step; private readonly GpuUInt128 step;
private readonly GpuUInt128 limit; private readonly GpuUInt128 limit;
PerfectNumbers.Core/HeuristicPrimeTester.cs
+47 -60
View File
@@ -86,15 +86,15 @@ public sealed class HeuristicPrimeTester
ct.ThrowIfCancellationRequested(); ct.ThrowIfCancellationRequested();
byte nMod10 = (byte)n.Mod10(); byte nMod10 = (byte)n.Mod10();
ulong sqrtLimit = ComputeHeuristicSqrt(n); ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
if (sqrtLimit < 3UL) if (maxDivisorSquare < 9UL)
{ {
return EvaluateWithOpenNumericFallback(n); return EvaluateWithOpenNumericFallback(n);
} }
bool includeGroupB = UseHeuristicGroupBTrialDivision; bool includeGroupB = UseHeuristicGroupBTrialDivision;
return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB);
} }
public bool IsPrimeCpu(ulong n) public bool IsPrimeCpu(ulong n)
@@ -105,33 +105,33 @@ public sealed class HeuristicPrimeTester
public bool IsPrimeGpu(ulong n) public bool IsPrimeGpu(ulong n)
{ {
byte nMod10 = (byte)n.Mod10(); byte nMod10 = (byte)n.Mod10();
ulong sqrtLimit = ComputeHeuristicSqrt(n); ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
return IsPrimeGpu(n, sqrtLimit, nMod10); return IsPrimeGpu(n, maxDivisorSquare, nMod10);
} }
public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) public bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10)
{ {
return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); return HeuristicIsPrimeGpuCore(n, maxDivisorSquare, nMod10);
} }
private bool HeuristicIsPrimeGpuCore(ulong n, ulong sqrtLimit, byte nMod10) private bool HeuristicIsPrimeGpuCore(ulong n, ulong maxDivisorSquare, byte nMod10)
{ {
// TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
if (sqrtLimit < 3UL) if (maxDivisorSquare < 9UL)
{ {
return EvaluateWithOpenNumericFallback(n); return EvaluateWithOpenNumericFallback(n);
} }
// if (!UseHeuristicGroupBTrialDivision) // if (!UseHeuristicGroupBTrialDivision)
// { // {
// return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); // return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB: false);
// } // }
bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10);
return !compositeDetected; return !compositeDetected;
} }
private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB) private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10, bool includeGroupB)
{ {
ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors; ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize);
@@ -145,7 +145,7 @@ public sealed class HeuristicPrimeTester
while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
{ {
ulong divisor = (ulong)groupADivisors[groupAIndex]; ulong divisor = (ulong)groupADivisors[groupAIndex];
if (divisor > sqrtLimit) if (!DivisorWithinSquareLimit(divisor, maxDivisorSquare))
{ {
groupAIndex = groupADivisors.Length; groupAIndex = groupADivisors.Length;
break; break;
@@ -165,7 +165,7 @@ public sealed class HeuristicPrimeTester
break; break;
} }
if ((ulong)groupADivisors[groupAIndex] > sqrtLimit) if (!DivisorWithinSquareLimit((ulong)groupADivisors[groupAIndex], maxDivisorSquare))
{ {
groupAIndex = groupADivisors.Length; groupAIndex = groupADivisors.Length;
break; break;
@@ -185,8 +185,8 @@ public sealed class HeuristicPrimeTester
? stackalloc int[endingOrder.Length] ? stackalloc int[endingOrder.Length]
: new int[endingOrder.Length]; : new int[endingOrder.Length];
bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; bool groupAHasMore = groupAIndex < groupADivisors.Length && DivisorWithinSquareLimit((ulong)groupADivisors[groupAIndex], maxDivisorSquare);
bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
while (groupAHasMore || groupBHasMore) while (groupAHasMore || groupBHasMore)
{ {
@@ -196,7 +196,7 @@ public sealed class HeuristicPrimeTester
while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
{ {
ulong divisor = (ulong)groupADivisors[groupAIndex]; ulong divisor = (ulong)groupADivisors[groupAIndex];
if (divisor > sqrtLimit) if (!DivisorWithinSquareLimit(divisor, maxDivisorSquare))
{ {
groupAIndex = groupADivisors.Length; groupAIndex = groupADivisors.Length;
groupAHasMore = false; groupAHasMore = false;
@@ -218,7 +218,7 @@ public sealed class HeuristicPrimeTester
} }
else else
{ {
groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; groupAHasMore = DivisorWithinSquareLimit((ulong)groupADivisors[groupAIndex], maxDivisorSquare);
if (!groupAHasMore) if (!groupAHasMore)
{ {
groupAIndex = groupADivisors.Length; groupAIndex = groupADivisors.Length;
@@ -231,7 +231,7 @@ public sealed class HeuristicPrimeTester
int processed = 0; int processed = 0;
while (processed < interleaveBatchSize) while (processed < interleaveBatchSize)
{ {
if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor))
{ {
groupBHasMore = false; groupBHasMore = false;
break; break;
@@ -247,7 +247,7 @@ public sealed class HeuristicPrimeTester
if (groupBHasMore) if (groupBHasMore)
{ {
groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
} }
} }
} }
@@ -258,7 +258,7 @@ public sealed class HeuristicPrimeTester
private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
{ {
int batchCapacity = HeuristicGpuDivisorBatchSize; int batchCapacity = HeuristicGpuDivisorBatchSize;
int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize; int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize;
@@ -296,8 +296,8 @@ public sealed class HeuristicPrimeTester
? stackalloc int[endingOrder.Length] ? stackalloc int[endingOrder.Length]
: new int[endingOrder.Length]; : new int[endingOrder.Length];
bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; bool groupAHasMore = groupAIndex < groupADivisors.Length && DivisorWithinSquareLimit((ulong)groupADivisors[groupAIndex], maxDivisorSquare);
bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, sqrtLimit); bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
while (!compositeDetected && (groupAHasMore || groupBHasMore)) while (!compositeDetected && (groupAHasMore || groupBHasMore))
{ {
@@ -307,7 +307,7 @@ public sealed class HeuristicPrimeTester
while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
{ {
ulong divisor = (ulong)groupADivisors[groupAIndex]; ulong divisor = (ulong)groupADivisors[groupAIndex];
if (divisor > sqrtLimit) if (!DivisorWithinSquareLimit(divisor, maxDivisorSquare))
{ {
groupAIndex = groupADivisors.Length; groupAIndex = groupADivisors.Length;
groupAHasMore = false; groupAHasMore = false;
@@ -338,7 +338,7 @@ public sealed class HeuristicPrimeTester
} }
else else
{ {
groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; groupAHasMore = DivisorWithinSquareLimit((ulong)groupADivisors[groupAIndex], maxDivisorSquare);
if (!groupAHasMore) if (!groupAHasMore)
{ {
groupAIndex = groupADivisors.Length; groupAIndex = groupADivisors.Length;
@@ -351,7 +351,7 @@ public sealed class HeuristicPrimeTester
int processed = 0; int processed = 0;
while (processed < interleaveBatchSize) while (processed < interleaveBatchSize)
{ {
if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor))
{ {
groupBHasMore = false; groupBHasMore = false;
break; break;
@@ -376,7 +376,7 @@ public sealed class HeuristicPrimeTester
if (groupBHasMore) if (groupBHasMore)
{ {
groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
} }
} }
} }
@@ -399,28 +399,15 @@ Cleanup:
{ {
return Prime.Numbers.IsPrime(n); return Prime.Numbers.IsPrime(n);
} }
internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n)
internal static ulong ComputeHeuristicSqrt(ulong n)
{ {
ulong sqrt = (ulong)Math.Sqrt(n); return n;
UInt128 square = (UInt128)sqrt * sqrt; }
while (square > n) [MethodImpl(MethodImplOptions.AggressiveInlining)]
{ private static bool DivisorWithinSquareLimit(ulong divisor, ulong maxDivisorSquare)
sqrt--; {
square = (UInt128)sqrt * sqrt; return divisor != 0UL && divisor <= maxDivisorSquare / divisor;
}
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 private static ReadOnlySpan<byte> GetGroupBEndingOrder(byte nMod10) => nMod10 switch
@@ -434,7 +421,7 @@ Cleanup:
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong sqrtLimit, out ulong divisor) private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong maxDivisorSquare, out ulong divisor)
{ {
ulong bestCandidate = ulong.MaxValue; ulong bestCandidate = ulong.MaxValue;
int bestEndingIndex = -1; int bestEndingIndex = -1;
@@ -449,7 +436,7 @@ Cleanup:
} }
ulong candidate = (ulong)divisors[index]; ulong candidate = (ulong)divisors[index];
if (candidate > sqrtLimit) if (!DivisorWithinSquareLimit(candidate, maxDivisorSquare))
{ {
indices[i] = divisors.Length; indices[i] = divisors.Length;
continue; continue;
@@ -474,7 +461,7 @@ Cleanup:
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool HasGroupBCandidates(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong sqrtLimit) private static bool HasGroupBCandidates(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong maxDivisorSquare)
{ {
for (int i = 0; i < endingOrder.Length; i++) for (int i = 0; i < endingOrder.Length; i++)
{ {
@@ -486,7 +473,7 @@ Cleanup:
} }
ulong candidate = (ulong)divisors[index]; ulong candidate = (ulong)divisors[index];
if (candidate > sqrtLimit) if (!DivisorWithinSquareLimit(candidate, maxDivisorSquare))
{ {
indices[i] = divisors.Length; indices[i] = divisors.Length;
continue; continue;
@@ -533,9 +520,9 @@ Cleanup:
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer) internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
{ {
return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer);
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -613,20 +600,20 @@ Cleanup:
internal ref struct HeuristicDivisorEnumerator internal ref struct HeuristicDivisorEnumerator
{ {
private readonly ulong sqrtLimit; private readonly ulong maxDivisorSquare;
private int groupAConstantIndex; private int groupAConstantIndex;
private ulong groupAWheelCandidate; private ulong groupAWheelCandidate;
private int groupAIncrementIndex; private int groupAIncrementIndex;
private bool groupAWheelActive; private bool groupAWheelActive;
private Span<HeuristicGroupBSequenceState> groupBStates; private Span<HeuristicGroupBSequenceState> groupBStates;
public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer) public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
{ {
this.sqrtLimit = sqrtLimit; this.maxDivisorSquare = maxDivisorSquare;
groupAConstantIndex = 0; groupAConstantIndex = 0;
groupAWheelCandidate = 23UL; groupAWheelCandidate = 23UL;
groupAIncrementIndex = 0; groupAIncrementIndex = 0;
groupAWheelActive = sqrtLimit >= 23UL; groupAWheelActive = DivisorWithinSquareLimit(23UL, maxDivisorSquare);
int count = InitializeGroupBStates(nMod10, groupBBuffer); int count = InitializeGroupBStates(nMod10, groupBBuffer);
groupBStates = count == 0 ? Span<HeuristicGroupBSequenceState>.Empty : groupBBuffer[..count]; groupBStates = count == 0 ? Span<HeuristicGroupBSequenceState>.Empty : groupBBuffer[..count];
@@ -643,7 +630,7 @@ Cleanup:
while (groupAConstantIndex < GroupAConstantDivisors.Length) while (groupAConstantIndex < GroupAConstantDivisors.Length)
{ {
ulong value = GroupAConstantDivisors[groupAConstantIndex++]; ulong value = GroupAConstantDivisors[groupAConstantIndex++];
if (value > sqrtLimit) if (!DivisorWithinSquareLimit(value, maxDivisorSquare))
{ {
groupAConstantIndex = GroupAConstantDivisors.Length; groupAConstantIndex = GroupAConstantDivisors.Length;
break; break;
@@ -657,7 +644,7 @@ Cleanup:
if (groupAWheelActive) if (groupAWheelActive)
{ {
ulong value = groupAWheelCandidate; ulong value = groupAWheelCandidate;
if (value <= sqrtLimit) if (DivisorWithinSquareLimit(value, maxDivisorSquare))
{ {
candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAWheel, 3, 0, (ushort)(value % Wheel210)); candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAWheel, 3, 0, (ushort)(value % Wheel210));
byte increment = GroupAIncrementPattern[groupAIncrementIndex]; byte increment = GroupAIncrementPattern[groupAIncrementIndex];
@@ -691,7 +678,7 @@ Cleanup:
continue; continue;
} }
if (value > sqrtLimit) if (!DivisorWithinSquareLimit(value, maxDivisorSquare))
{ {
state.Candidate = ulong.MaxValue; state.Candidate = ulong.MaxValue;
continue; continue;
docs/HeuristicIsPrimePlan.md
+1 -1
View File
@@ -131,7 +131,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels,
* [done] The CPU heuristic now operates as a pure trial-division sweep without maintaining per-call summaries, relying on `PrepareHeuristicDivisor` only when downstream consumers request Montgomery data. * [done] The CPU heuristic now operates as a pure trial-division sweep without maintaining per-call summaries, relying on `PrepareHeuristicDivisor` only when downstream consumers request Montgomery data.
* [done] Introduced the temporary `UseHeuristicGroupBTrialDivision` gate so current builds execute Group A locally and then fall back to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, keeping the full Group B implementation available for reactivation. * [done] Introduced the temporary `UseHeuristicGroupBTrialDivision` gate so current builds execute Group A locally and then fall back to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, keeping the full Group B implementation available for reactivation.
7. [done] **Implement `IsPrimeGpu`** 7. [done] **Implement `IsPrimeGpu`**
* [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `sqrtLimit`/`nMod10` inputs, and falling back to CPU when GPUs are disabled. * [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `maxDivisorSquare`/`nMod10` inputs, and falling back to CPU when GPUs are disabled.
* [done] Added `PrimeTesterKernels.HeuristicTrialDivisionKernel` so GPU batches respect Group A/B ordering while emitting hit flags for early exits. * [done] Added `PrimeTesterKernels.HeuristicTrialDivisionKernel` so GPU batches respect Group A/B ordering while emitting hit flags for early exits.
* [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic. * [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic.
* [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement. * [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement.