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

Use heuristic sieves for prime trial division

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-10-26 03:41:04 +01:00
1 parent 15c229665d
commit 1f7680a465
4 files changed
+235 -211
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/HeuristicPrimeSieves.cs
+79 -163
View File
@@ -15,195 +15,111 @@ public static class HeuristicPrimeSieves
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 readonly int[] GroupADivisorsStorage;
private static readonly ulong[] GroupADivisorSquaresStorage;
private static int[]? s_groupBDivisorsEnding1;
private static ulong[]? s_groupBDivisorSquaresEnding1;
private static readonly int[] GroupBDivisorsEnding1Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding1Storage;
private static int[]? s_groupBDivisorsEnding3;
private static ulong[]? s_groupBDivisorSquaresEnding3;
private static readonly int[] GroupBDivisorsEnding3Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding3Storage;
private static int[]? s_groupBDivisorsEnding7;
private static ulong[]? s_groupBDivisorSquaresEnding7;
private static readonly int[] GroupBDivisorsEnding7Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding7Storage;
private static int[]? s_groupBDivisorsEnding9;
private static ulong[]? s_groupBDivisorSquaresEnding9;
private static readonly int[] GroupBDivisorsEnding9Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding9Storage;
private static readonly object InitLock = new();
public static ReadOnlySpan<int> GroupADivisors
static HeuristicPrimeSieves()
{
get
{
EnsureGroupAInitialized();
return s_groupADivisors!;
}
Console.WriteLine("[HeuristicPrimeSieves] Initializing heuristic divisor sieves...");
(GroupADivisorsStorage, GroupADivisorSquaresStorage) = BuildGroupADivisors();
(GroupBDivisorsEnding1Storage, GroupBDivisorSquaresEnding1Storage) = BuildGroupBResidueSieve(GroupBResiduesEnding1);
(GroupBDivisorsEnding3Storage, GroupBDivisorSquaresEnding3Storage) = BuildGroupBResidueSieve(GroupBResiduesEnding3);
(GroupBDivisorsEnding7Storage, GroupBDivisorSquaresEnding7Storage) = BuildGroupBResidueSieve(GroupBResiduesEnding7);
(GroupBDivisorsEnding9Storage, GroupBDivisorSquaresEnding9Storage) = BuildGroupBResidueSieve(GroupBResiduesEnding9);
}
public static ReadOnlySpan<ulong> GroupADivisorSquares
public static void EnsureInitialized()
{
get
{
EnsureGroupAInitialized();
return s_groupADivisorSquares!;
}
// Intentionally left blank. Accessing this method forces the static constructor to run.
}
public static ReadOnlySpan<int> GroupBDivisorsEnding1
public static ReadOnlySpan<int> GroupADivisors => GroupADivisorsStorage;
public static ReadOnlySpan<ulong> GroupADivisorSquares => GroupADivisorSquaresStorage;
public static ReadOnlySpan<int> GroupBDivisorsEnding1 => GroupBDivisorsEnding1Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding1 => GroupBDivisorSquaresEnding1Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding3 => GroupBDivisorsEnding3Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding3 => GroupBDivisorSquaresEnding3Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding7 => GroupBDivisorsEnding7Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding7 => GroupBDivisorSquaresEnding7Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding9 => GroupBDivisorsEnding9Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding9 => GroupBDivisorSquaresEnding9Storage;
private static (int[] divisors, ulong[] squares) BuildGroupADivisors()
{
get
var divisors = new int[SieveLength];
var squares = new ulong[SieveLength];
int index = 0;
for (; index < GroupAConstants.Length && index < divisors.Length; index++)
{
EnsureGroupBInitialized(ref s_groupBDivisorsEnding1, ref s_groupBDivisorSquaresEnding1, GroupBResiduesEnding1);
return s_groupBDivisorsEnding1!;
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++;
}
return (divisors, squares);
}
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding1
private static (int[] divisors, ulong[] squares) BuildGroupBResidueSieve(ReadOnlySpan<ushort> residues)
{
get
{
EnsureGroupBInitialized(ref s_groupBDivisorsEnding1, ref s_groupBDivisorSquaresEnding1, GroupBResiduesEnding1);
return s_groupBDivisorSquaresEnding1!;
}
}
var divisors = new int[SieveLength];
var squares = new ulong[SieveLength];
public static ReadOnlySpan<int> GroupBDivisorsEnding3
{
get
{
EnsureGroupBInitialized(ref s_groupBDivisorsEnding3, ref s_groupBDivisorSquaresEnding3, GroupBResiduesEnding3);
return s_groupBDivisorsEnding3!;
}
}
int index = 0;
ulong baseValue = 0UL;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding3
{
get
while (index < divisors.Length)
{
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)
for (int i = 0; i < residues.Length && index < divisors.Length; i++)
{
return;
}
ulong candidate = baseValue + residues[i];
if (candidate <= 13UL)
{
continue;
}
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;
divisors[index] = (int)candidate;
squares[index] = candidate * candidate;
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;
baseValue += (ulong)Wheel210;
}
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;
}
return (divisors, squares);
}
}
PerfectNumbers.Core/PrimeTester.cs
+148 -40
View File
@@ -14,6 +14,11 @@ namespace PerfectNumbers.Core;
public sealed class PrimeTester(bool useInternal = false)
{
static PrimeTester()
{
HeuristicPrimeSieves.EnsureInitialized();
}
// 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];
@@ -197,20 +202,13 @@ public sealed class PrimeTester(bool useInternal = false)
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))
ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
for (int i = 0; i < groupADivisors.Length; i++)
{
if (!includeGroupB && candidate.Group == HeuristicDivisorGroup.GroupB)
ulong divisor = (ulong)groupADivisors[i];
if (divisor > sqrtLimit)
{
return EvaluateWithOpenNumericFallback(n);
}
ulong divisor = candidate.Value;
if (divisor <= 1UL)
{
continue;
break;
}
if (n % divisor == 0UL)
@@ -219,20 +217,71 @@ public sealed class PrimeTester(bool useInternal = false)
}
}
return includeGroupB ? true : EvaluateWithOpenNumericFallback(n);
if (!includeGroupB)
{
return EvaluateWithOpenNumericFallback(n);
}
ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
if (endingOrder.IsEmpty)
{
return EvaluateWithOpenNumericFallback(n);
}
Span<int> indices = endingOrder.Length <= 8
? stackalloc int[endingOrder.Length]
: new int[endingOrder.Length];
while (true)
{
ulong bestCandidate = ulong.MaxValue;
int bestEndingIndex = -1;
for (int i = 0; i < endingOrder.Length; i++)
{
ReadOnlySpan<int> divisors = GetGroupBDivisors(endingOrder[i]);
int index = indices[i];
if ((uint)index >= (uint)divisors.Length)
{
continue;
}
ulong candidate = (ulong)divisors[index];
if (candidate > sqrtLimit)
{
indices[i] = divisors.Length;
continue;
}
if (bestEndingIndex == -1 || candidate < bestCandidate)
{
bestCandidate = candidate;
bestEndingIndex = i;
}
}
if (bestEndingIndex == -1)
{
break;
}
if (n % bestCandidate == 0UL)
{
return false;
}
indices[bestEndingIndex]++;
}
return true;
}
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;
@@ -243,7 +292,6 @@ public sealed class PrimeTester(bool useInternal = false)
try
{
candidateArray = candidatePool.Rent(batchCapacity);
divisorArray = divisorPool.Rent(batchCapacity);
hitFlags = hitPool.Rent(batchCapacity);
@@ -259,6 +307,11 @@ public sealed class PrimeTester(bool useInternal = false)
bool ProcessBatch(int length)
{
if (length <= 0)
{
return false;
}
scratch.Input.View.CopyFromCPU(ref divisorArray![0], length);
state.HeuristicTrialDivisionKernel(length, scratch.Input.View, n, scratch.Output.View);
accelerator.Synchronize();
@@ -271,7 +324,7 @@ public sealed class PrimeTester(bool useInternal = false)
continue;
}
ulong divisor = candidateArray![i].Value;
ulong divisor = divisorArray![i];
if (divisor > 1UL && n % divisor == 0UL)
{
return true;
@@ -281,36 +334,87 @@ public sealed class PrimeTester(bool useInternal = false)
return false;
}
while (enumerator.TryGetNext(out HeuristicDivisorCandidate candidate))
bool QueueDivisor(ulong divisor)
{
ulong divisor = candidate.Value;
if (divisor <= 1UL)
{
continue;
}
candidateArray[count] = candidate;
divisorArray[count] = divisor;
divisorArray![count] = divisor;
count++;
if (count == batchCapacity)
{
if (ProcessBatch(count))
{
compositeDetected = true;
break;
return true;
}
count = 0;
}
return false;
}
ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
for (int i = 0; i < groupADivisors.Length && !compositeDetected; i++)
{
ulong divisor = (ulong)groupADivisors[i];
if (divisor > sqrtLimit)
{
break;
}
compositeDetected = QueueDivisor(divisor);
}
if (!compositeDetected)
{
ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
if (!endingOrder.IsEmpty)
{
Span<int> indices = endingOrder.Length <= 8
? stackalloc int[endingOrder.Length]
: new int[endingOrder.Length];
while (!compositeDetected)
{
ulong bestCandidate = ulong.MaxValue;
int bestEndingIndex = -1;
for (int i = 0; i < endingOrder.Length; i++)
{
ReadOnlySpan<int> divisors = GetGroupBDivisors(endingOrder[i]);
int index = indices[i];
if ((uint)index >= (uint)divisors.Length)
{
continue;
}
ulong candidate = (ulong)divisors[index];
if (candidate > sqrtLimit)
{
indices[i] = divisors.Length;
continue;
}
if (bestEndingIndex == -1 || candidate < bestCandidate)
{
bestCandidate = candidate;
bestEndingIndex = i;
}
}
if (bestEndingIndex == -1)
{
break;
}
compositeDetected = QueueDivisor(bestCandidate);
indices[bestEndingIndex]++;
}
}
}
if (!compositeDetected && count > 0)
{
if (ProcessBatch(count))
{
compositeDetected = true;
}
compositeDetected = ProcessBatch(count);
}
}
finally
@@ -331,11 +435,6 @@ public sealed class PrimeTester(bool useInternal = false)
divisorPool.Return(divisorArray);
}
if (candidateArray is not null)
{
candidatePool.Return(candidateArray);
}
gpu.Dispose();
limiter.Dispose();
}
@@ -380,6 +479,15 @@ public sealed class PrimeTester(bool useInternal = false)
_ => ReadOnlySpan<byte>.Empty,
};
private static ReadOnlySpan<int> GetGroupBDivisors(byte ending) => ending switch
{
1 => HeuristicPrimeSieves.GroupBDivisorsEnding1,
3 => HeuristicPrimeSieves.GroupBDivisorsEnding3,
7 => HeuristicPrimeSieves.GroupBDivisorsEnding7,
9 => HeuristicPrimeSieves.GroupBDivisorsEnding9,
_ => ReadOnlySpan<int>.Empty,
};
private static ReadOnlySpan<ushort> GetWheel210ResiduesForEnding(byte ending) => ending switch
{
1 => Wheel210ResiduesEnding1,