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

Merge pull request #222 from mlemanczyk/codex/refactor-primetester.isprimegpu-for-branch-less-execution

Browse Source 
Clarify GPU sieve invariants and trim defensive guards
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-10-22 05:24:57 +02:00
2 parents ab4388dfba f2070a0f64
commit 6c3efb90c9
4 files changed
+108 -103
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs
+51 -64
View File
@@ -10,52 +10,38 @@ internal static class PrimeTesterKernels
public static void SmallPrimeSieveKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<uint> smallPrimes, ArrayView<byte> results)
{
ulong n = numbers[index];
if (n <= 3UL)
// EvenPerfectBitScanner filters candidates so production GPU launches only see:
// - n >= 31 (the scanner starts at 136,279,841 and test mode seeds 31)
// - n odd (candidate generators and ModResidueTracker eliminate even exponents)
// - n not divisible by 5 (residue prefiltering and AddPrimes transform cover these composites)
// Tests should honor the same constraints so the host wrapper can remain branchless.
ulong log2Estimate = 63UL - (ulong)XMath.LeadingZeroCount(n | 1UL);
ulong gcd = BinaryGcdGpu(n, log2Estimate);
bool heuristicActive = n.Mod10() == 1UL;
bool gcdReject = heuristicActive & (gcd != 1UL);
bool composite = gcdReject;
bool continueLoop = !gcdReject;
long length = smallPrimes.Length;
for (int i = 0; i < length; i++)
{
results[index] = (byte)(n >= 2UL ? 1 : 0);
return;
ulong prime = smallPrimes[i];
ulong primeSquare = (ulong)prime * prime;
bool withinRange = primeSquare <= n;
ulong remainder = n % prime;
bool divides = remainder == 0UL;
bool markComposite = continueLoop & withinRange & divides;
composite |= markComposite;
continueLoop &= withinRange & !divides;
}
if ((n & 1UL) == 0UL || (n > 5UL && (n % 5UL) == 0UL))
{
// TODO: Replace the `% 5` branch with the GPU Mod5 helper once the sieve kernel can
// reuse the benchmarked bitmask reduction instead of modulo on every candidate.
results[index] = 0;
return;
}
if (n.Mod10() == 1UL)
{
// Early reject special GCD heuristic with floor(log2 n)
ulong m = 63UL - (ulong)XMath.LeadingZeroCount(n);
if (BinaryGcdGpu(n, m) != 1UL)
{
results[index] = 0;
return;
}
}
long len = smallPrimes.Length;
for (int i = 0; i < len; i++)
{
ulong p = smallPrimes[i];
if (p * p > n)
{
break;
}
if ((n % p) == 0UL)
{
// TODO: Route this modulo through the shared divisor-cycle cache once exposed to GPU kernels
// so batched sieves avoid per-prime `%` operations that profiling showed expensive.
results[index] = 0;
return;
}
}
results[index] = 1;
results[index] = composite ? (byte)0 : (byte)1;
}
public static void SharesFactorKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<byte> results)
{
ulong n = numbers[index];
@@ -69,31 +55,32 @@ internal static class PrimeTesterKernels
// TODO: Replace this inline GPU binary GCD with the kernel extracted from
// GpuUInt128BinaryGcdBenchmarks via GpuKernelPool so device callers reuse the
// fully unrolled ladder instead of this branchy fallback.
if (u == 0UL)
bool uIsZero = u == 0UL;
bool vIsZero = v == 0UL;
bool eitherZero = uIsZero | vIsZero;
ulong zeroResult = uIsZero ? v : u;
ulong combined = u | v | 1UL;
int shift = XMath.TrailingZeroCount(combined);
ulong normalizedU = u >> XMath.TrailingZeroCount(u | 1UL);
ulong normalizedV = v >> XMath.TrailingZeroCount(v | 1UL);
ulong currentU = normalizedU;
ulong currentV = normalizedV;
bool loopCondition = currentV != 0UL;
while (loopCondition)
{
return v;
currentV >>= XMath.TrailingZeroCount(currentV);
ulong minValue = XMath.Min(currentU, currentV);
ulong maxValue = XMath.Max(currentU, currentV);
currentU = minValue;
currentV = maxValue - minValue;
loopCondition = currentV != 0UL;
}
if (v == 0UL)
{
return u;
}
int shift = XMath.TrailingZeroCount(u | v);
u >>= XMath.TrailingZeroCount(u);
do
{
v >>= XMath.TrailingZeroCount(v);
if (u > v)
{
(u, v) = (v, u);
}
v -= u;
}
while (v != 0UL);
return u << shift;
ulong gcd = currentU << shift;
return eitherZero ? zeroResult : gcd;
}
}
PerfectNumbers.Core/PrimeTester.cs
+51 -35
View File
@@ -16,61 +16,77 @@ public sealed class PrimeTester(bool useInternal = false)
// indirection shows up when we sieve millions of candidates per second.
IsPrimeInternal(n, ct);
public static bool IsPrimeGpu(ulong n) => new PrimeTester().IsPrimeGpu(n, CancellationToken.None);
public static bool IsPrimeGpu(ulong n) => new PrimeTester().IsPrimeGpu(n, CancellationToken.None);
// Optional GPU-assisted primality: batched small-prime sieve on device.
public bool IsPrimeGpu(ulong n, CancellationToken ct)
// Optional GPU-assisted primality: batched small-prime sieve on device.
public bool IsPrimeGpu(ulong n, CancellationToken ct)
{
if (GpuContextPool.ForceCpu)
{
return IsPrimeInternal(n, 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.
// if (n <= 3UL)
// {
// return n >= 2UL;
// }
// Candidate generation (Add/AddPrimes transforms plus ModResidueTracker) already strips even values and multiples of five.
// The downstream small-factor sweep and Pollard routines factor divisors q = 2kp + 1 with k >= 1 without feeding new
// exponents back into this path, so the GPU sees the same odd sequence that already passed the CPU residues.
// 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 scanner pipeline.
// if (n < 31UL)
// {
// return IsPrimeInternal(n, ct);
// }
bool forceCpu = GpuContextPool.ForceCpu;
Span<ulong> one = stackalloc ulong[1];
Span<byte> outFlags = stackalloc byte[1];
// TODO: Inline the single-value GPU sieve fast path from GpuModularArithmeticBenchmarks so this wrapper
// can skip stackalloc buffers and reuse the pinned upload span the benchmark identified as fastest.
one[0] = n;
IsPrimeBatchGpu(one, outFlags);
outFlags[0] = 0;
if (outFlags[0] != 0)
if (!forceCpu)
{
return true;
IsPrimeBatchGpu(one, outFlags);
}
// Defensive fallback: GPU sieve produced a composite verdict.
// Confirm with CPU logic to prevent false negatives observed on
// some accelerators.
return IsPrimeInternal(n, ct);
bool belowGpuRange = n < 31UL;
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;
}
// 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)
{
bool result = true;
if (n <= 3UL)
// 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
// this path for arbitrary residues and cofactors, so keep the boolean guards even though the scanner never trips them.
bool isTwo = n == 2UL;
bool isOdd = (n & 1UL) != 0UL;
// TODO: Replace this modulo check with ULongExtensions.Mod5 so the CPU hot path reuses the benchmarked helper instead of `%`.
ulong mod5 = n % 5UL;
bool divisibleByFive = n > 5UL && mod5 == 0UL;
bool result = n >= 2UL && (isTwo || isOdd) && !divisibleByFive;
bool requiresTrialDivision = result && n >= 7UL && !isTwo;
if (requiresTrialDivision)
{
result = n >= 2UL;
}
else if ((n & 1UL) == 0)
{
result = false;
}
else if (n > 5UL && (n % 5UL) == 0UL)
{
// TODO: Replace this modulo check with ULongExtensions.Mod5 so the CPU hot path
// reuses the benchmarked helper instead of `%` when sieving large prime ranges.
result = false;
}
else
{
if (n.Mod10() == 1UL && SharesFactorWithMaxExponent(n))
{
result = false;
}
else
bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n);
result &= !sharesMaxExponentFactor;
if (result)
{
var smallPrimeDivisorsLength = PrimesGenerator.SmallPrimes.Length;
uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes;