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

Clarify cycle stepper reset inlining

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-10-02 13:30:56 +02:00
1 parent 8b467c27f3
commit b8aa845c5c
53 files changed
+1736 -1131
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/Gpu/GpuConstants.cs
+8 -5
View File
@@ -1,8 +1,11 @@
namespace PerfectNumbers.Core.Gpu
{
public static class GpuConstants
{
public const int ScanBatchSize = 2_097_152;
public const int GpuCycleStepsPerInvocation = 256;
public static class GpuConstants
{
public const int ScanBatchSize = 2_097_152;
public const int GpuCycleStepsPerInvocation = 256;
// TODO: Auto-tune these GPU constants per accelerator so batch sizing aligns with the fastest kernel
// configuration reported by the Pow2Mod and divisor-cycle benchmarks.
}
}
}
PerfectNumbers.Core/Gpu/GpuContextPool.cs
+6 -3
View File
@@ -6,9 +6,11 @@ namespace PerfectNumbers.Core.Gpu;
public static class GpuContextPool
{
private static readonly bool PoolingEnabled = true;
// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
public static bool ForceCpu { get; set; } = false;
private static readonly bool PoolingEnabled = true;
// TODO: Introduce accelerator-specific warmup so pooled contexts precompile the ProcessEightBitWindows kernels and load
// divisor-cycle data before the first scan begins.
// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
public static bool ForceCpu { get; set; } = false;
internal sealed class PooledContext
{
@@ -156,3 +158,4 @@ public static class GpuContextPool
}
}
}
PerfectNumbers.Core/Gpu/GpuKernelPool.cs
+4 -4
View File
@@ -694,7 +694,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
public class GpuKernelPool
{
private static readonly ConcurrentDictionary<Accelerator, KernelContainer> KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, KernelContainer> KernelCache = new(); // TODO: Replace this concurrent map with a simple accelerator-indexed lookup once kernel launchers are prewarmed during startup so we can drop the thread-safe wrapper entirely.
private static KernelContainer GetKernels(Accelerator accelerator)
{
@@ -712,14 +712,14 @@ public class GpuKernelPool
}
// Ensure single upload per accelerator even if multiple threads race here.
lock (kernels)
lock (kernels) // TODO: Remove this lock by pre-uploading the immutable small-cycle snapshot during initialization; once no mutation happens at runtime, the pool must expose a simple reference without synchronization.
{
if (kernels.SmallCycles is { } existing)
{
return existing.View;
}
var host = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
var host = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); // TODO: Preload this device buffer during startup and keep it immutable so we can delete the lock above in favor of the preloaded snapshot.
var device = accelerator.Allocate1D<ulong>(host.Length);
device.View.CopyFromCPU(host);
kernels.SmallCycles = device;
@@ -738,7 +738,7 @@ public class GpuKernelPool
return new ResiduePrimeViews(lastOne.View, lastSeven.View, lastOnePow2.View, lastSevenPow2.View);
}
lock (kernels)
lock (kernels) // TODO: Inline these small-prime uploads into startup initialization alongside the small-cycle snapshot so we can drop runtime locking and keep the GPU pool free of synchronization.
{
if (kernels.SmallPrimesLastOne is { } existingLastOne &&
kernels.SmallPrimesLastSeven is { } existingLastSeven &&
PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs
+4
View File
@@ -8,6 +8,8 @@ public static class GpuPrimeWorkLimiter
public static IDisposable Acquire()
{
// TODO: Inline the pooled limiter guard from the GpuLimiterThroughputBenchmarks so prime-sieve
// GPU jobs stop allocating Releaser instances on every acquisition.
var sem = _semaphore;
sem.Wait();
return new Releaser(sem);
@@ -27,6 +29,8 @@ public static class GpuPrimeWorkLimiter
return;
}
// TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using
// the same semaphore pool without reallocating per adjustment.
_semaphore = new SemaphoreSlim(value, value);
_currentLimit = value;
}
PerfectNumbers.Core/Gpu/GpuUInt128.cs
+29 -16
View File
@@ -156,14 +156,14 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool operator >=(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) >= 0;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 operator ++(GpuUInt128 value)
{
value.Add(1UL);
return value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 operator ++(GpuUInt128 value)
{
value.Add(1UL);
return value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly int CompareTo(GpuUInt128 other)
{
if (High < other.High)
@@ -536,6 +536,9 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 BinaryGcd(GpuUInt128 u, GpuUInt128 v)
{
// TODO: Replace this scalar binary GCD with the branchless reduction from
// GpuUInt128BinaryGcdBenchmarks so CPU fallbacks stay aligned with the GPU kernel
// performance when resolving large divisor residues.
if (u.IsZero)
{
return v;
@@ -1044,11 +1047,15 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void ModPow(GpuUInt128 exponent, GpuUInt128 modulus)
{
// TODO: This should operate on the instance itself, not on a copy. Avoid creating new instances anywhere.
// TODO: Swap this copy-heavy path for the pooled base/exponent ladder used in
// GpuUInt128MulModBenchmarks so the GPU shim keeps reusing buffers instead of
// allocating temporary structs during residue scans.
GpuUInt128 result = new(1UL);
GpuUInt128 baseValue = this;
// TODO: Can we modify this loop to process multiple bits at a time? E.g. 64-bit chunks.
// TODO: Replace the single-bit square-and-multiply loop with the 64-bit windowed
// ladder measured fastest in GpuUInt128MulModBenchmarks to align with the GPU kernel
// implementation once the shared helper lands.
while (!exponent.IsZero)
{
if ((exponent.Low & 1UL) != 0UL)
@@ -1066,11 +1073,13 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void ModPow(ulong exponent, GpuUInt128 modulus)
{
// TODO: This should operate on the instance itself, not on a copy. Avoid creating new instances anywhere.
// TODO: Share the pooled ladder state from GpuUInt128MulModBenchmarks here as well so
// Lucas–Lehmer batches avoid constructing throwaway temporaries on every exponent.
GpuUInt128 result = new(1UL);
GpuUInt128 baseValue = this;
// TODO: Can we modify this loop to process multiple bits at a time? E.g. 64-bit chunks.
// TODO: Upgrade this loop to the same 64-bit windowed ladder proven fastest in
// GpuUInt128MulModBenchmarks so the scalar helper matches the GPU-optimized path.
while (exponent != 0UL)
{
if ((exponent & 1UL) != 0UL)
@@ -1086,14 +1095,16 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
Low = result.Low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void ModInv(ulong modulus)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void ModInv(ulong modulus)
{
// TODO: We shouldn't need a new instance here to represent modulus.
// TODO: Replace this Fermat inversion with the Montgomery ladder highlighted in
// GpuUInt128Montgomery64Benchmarks so we avoid instantiating a temporary modulus and
// reuse the pooled reduction constants.
ModPow(modulus - 2UL, new GpuUInt128(modulus));
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void ModInv(GpuUInt128 modulus)
{
if (modulus.High == 0UL)
@@ -1114,7 +1125,9 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void MultiplyFull(GpuUInt128 left, GpuUInt128 right, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
{
// TODO: Reuse variables to reduce register pressure.
// TODO: Reuse variables to reduce register pressure following the fused-limb layout
// from GpuUInt128MulModByLimbBenchmarks so the multiply helper matches the fastest
// GPU-compatible scalar routine.
var (h0, l0) = Mul64(left.Low, right.Low);
var (h1, l1) = Mul64(left.Low, right.High);
var (h2, l2) = Mul64(left.High, right.Low);
PerfectNumbers.Core/Gpu/GpuWorkLimiter.cs
+5
View File
@@ -8,6 +8,9 @@ public static class GpuWorkLimiter
public static IDisposable Acquire()
{
// TODO: Replace the per-call Releaser allocation with the pooled struct-based guard from the
// GpuLimiterThroughputBenchmarks so limiter acquisition does not allocate when we enter the
// GPU scanning hot path.
var sem = _semaphore;
sem.Wait();
return new Releaser(sem);
@@ -30,6 +33,8 @@ public static class GpuWorkLimiter
// Swap to a new semaphore for the new limit. We intentionally
// do not dispose the old semaphore to avoid releasing handles
// that may still be in use by existing Releaser instances.
// TODO: Consolidate with GpuPrimeWorkLimiter so both limiters share a pooled SemaphoreSlim and avoid rebuilding the
// limiter state whenever limits are adjusted from the CLI.
_semaphore = new SemaphoreSlim(value, value);
_currentLimit = value;
}
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+15
View File
@@ -369,6 +369,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
DivisorCycleCache divisorCyclesCache = DivisorCycleCache.Shared;
DivisorCycleCache.CycleBlock cycleLease = divisorCyclesCache.Acquire(currentDivisor);
// TODO: Replace the block lease with a direct single-cycle lookup so the GPU path
// follows the benchmarked single-block policy and computes the missing cycle inline on the
// configured device without mutating or swapping cache blocks after startup.
cycleLeaseEnd = cycleLease.End;
try
@@ -389,6 +392,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
{
cycleLease.Dispose();
cycleLease = divisorCyclesCache.Acquire(currentDivisor);
// TODO: Drop this block hand-off once the cache exposes a direct
// GetCycle helper that materializes only the requested divisor length
// without touching additional blocks or scheduling background work.
cycleLeaseEnd = cycleLease.End;
}
@@ -414,6 +420,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
exponent = 0UL;
}
}
else
{
// TODO: Trigger an immediate GPU computation for the single missing cycle without storing
// the result back or requesting extra cache blocks so very large divisors stay within the
// single-block constraint while still reusing the snapshot for other divisors.
}
if (includeDivisor)
{
@@ -617,6 +629,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
if (divisorCycle == 0UL)
{
// TODO: Replace this guard with an on-demand cycle computation that runs on the device selected
// for the session, keeps the result transient instead of caching it, and continues operating with the
// single shared block loaded from disk.
throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}.");
}
PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs
+4 -1
View File
@@ -110,7 +110,10 @@ public sealed class MersenneNumberDivisorGpuTester
continue;
}
UInt128 cycle128 = MersenneDivisorCycles.GetCycle(d);
UInt128 cycle128 = MersenneDivisorCycles.GetCycle(d);
// TODO: When this lookup falls outside the preloaded snapshot, invoke the single-cycle GPU helper
// instead of allowing GetCycle to allocate an entire block; the divisor-cycle benchmarks showed the
// one-off computation keeping memory flat while still delivering the optimized stepping data.
// TODO: Feed this check through the ProcessEightBitWindows-based residue tracker so we avoid the slow
// UInt128 modulo on the CPU fallback highlighted by the Pow2Montgomery benchmarks.
if ((UInt128)p % cycle128 != UInt128.Zero)
PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs
+89 -88
View File
@@ -6,98 +6,99 @@ namespace PerfectNumbers.Core.Gpu;
public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool useGpuOrder)
{
private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useGpuOrder = useGpuOrder;
private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useGpuOrder = useGpuOrder;
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope();
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy
UInt128 kStart = 1UL;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope();
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
int batchSize = GpuConstants.ScanBatchSize;
UInt128 kStart = 1UL;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
byte last = lastIsSeven ? (byte)1 : (byte)0;
var pow2Kernel = gpuLease.Pow2ModKernel; // TODO: Route this binding to the ProcessEightBitWindows kernel once the scalar helper lands so GPU incremental scans inherit the windowed speedup from GpuPow2ModBenchmarks.
var incKernel = gpuLease.IncrementalKernel;
exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
GpuUInt128 twoPGpu = (GpuUInt128)twoP;
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ResiduePrimeViews primeViews = default;
if (_kernelType == GpuKernelType.Pow2Mod)
{
primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
}
var pow2Kernel = gpuLease.Pow2ModKernel; // TODO: Route this binding to the ProcessEightBitWindows kernel once the scalar helper lands so GPU incremental scans inherit the windowed speedup from GpuPow2ModBenchmarks.
var incKernel = gpuLease.IncrementalKernel;
exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
GpuUInt128 twoPGpu = (GpuUInt128)twoP;
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ResiduePrimeViews primeViews = default;
if (_kernelType == GpuKernelType.Pow2Mod)
{
primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
}
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
// Avoid giant stack allocations that can trigger StackOverflow when batchSize is large.
// Rent a reusable array from the shared pool instead.
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
UInt128 remaining;
int currentSize;
int i;
UInt128 q = UInt128.Zero;
try
{
while (kStart <= maxK && Volatile.Read(ref isPrime))
{
remaining = maxK - kStart + UInt128.One;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize);
// Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + UInt128.One;
if (_kernelType == GpuKernelType.Pow2Mod)
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
pow2Kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul,
kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven,
primeViews.LastOnePow2, primeViews.LastSevenPow2);
}
else
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
incKernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul,
q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView);
}
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
UInt128 remaining;
int currentSize;
int i;
UInt128 q = UInt128.Zero;
stream.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
if (_kernelType == GpuKernelType.Pow2Mod)
{
for (i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
}
else
{
q = twoP * kStart + UInt128.One;
for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
{
if (orders[i] != 0UL && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
}
}
}
try
{
while (kStart <= maxK && Volatile.Read(ref isPrime))
{
remaining = maxK - kStart + UInt128.One;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize);
UInt128 q0 = twoP * kStart + UInt128.One;
// TODO: Swap this multiply-and-add initialization for the divisor-cycle stepping helper so the GPU incremental path reuses cached offsets instead of recomputing twoP multiples every batch.
kStart += (UInt128)currentSize;
}
}
finally
{
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
execution.Dispose();
gpuLease.Dispose();
}
if (_kernelType == GpuKernelType.Pow2Mod)
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
pow2Kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul,
kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven,
primeViews.LastOnePow2, primeViews.LastSevenPow2);
}
else
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
incKernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul,
q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView);
}
}
stream.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
if (_kernelType == GpuKernelType.Pow2Mod)
{
for (i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
}
else
{
q = twoP * kStart + UInt128.One;
// TODO: Replace this per-lane multiplication and the subsequent q += twoP stride with the shared divisor-cycle ladder so each iteration advances by cached cycle increments and automatically triggers on-demand GPU computations when the snapshot lacks the requested cycle.
for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
{
if (orders[i] != 0UL && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
}
}
}
kStart += (UInt128)currentSize;
}
}
finally
{
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
execution.Dispose();
gpuLease.Dispose();
}
}
}
PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs
+17
View File
@@ -43,6 +43,9 @@ public class MersenneNumberLucasLehmerGpuTester
{
// Default to global kernel preference for backward compatibility.
bool runOnGpu = !GpuContextPool.ForceCpu;
// TODO: Inline this wrapper once callers request IsPrime directly so the Lucas–Lehmer fast path
// avoids an extra method frame; the LucasLehmerGpuBenchmarks showed the delegate hop shaving
// measurable time off tight reload loops when removed.
return IsPrime(exponent, runOnGpu);
}
@@ -120,6 +123,8 @@ public class MersenneNumberLucasLehmerGpuTester
GpuUInt128[] modulusArray = ArrayPool<GpuUInt128>.Shared.Rent(count);
for (int i = 0; i < count; i++)
{
// TODO: Replace this per-exponent shift with a small shared table (one entry per supported exponent < 128)
// so Lucas–Lehmer batch runs reuse cached GpuUInt128 moduli instead of rebuilding them for every request.
modulusArray[i] = new GpuUInt128(((UInt128)1 << (int)exponents[i]) - 1UL);
}
@@ -292,11 +297,15 @@ public class MersenneNumberLucasLehmerGpuTester
continue;
}
// TODO: Replace int.TryParse with the span-based Utf8Parser helper when loading cached parameters so we avoid
// culture-aware parsing in this hot startup loop.
if (!int.TryParse(parts[0], out var storedLength) || storedLength != length)
{
continue;
}
// TODO: Switch these ulong.Parse calls to the Utf8Parser-based fast-path once we expose a zero-allocation reader for
// persisted kernel parameters.
ulong modHigh = ulong.Parse(parts[1]);
ulong modLow = ulong.Parse(parts[2]);
ulong rootHigh = ulong.Parse(parts[3]);
@@ -315,6 +324,9 @@ public class MersenneNumberLucasLehmerGpuTester
{
lock (ParameterFileLock)
{
// TODO: Replace StreamWriter with the pooled TextFileWriter pipeline so persisting NTT
// parameters reuses the zero-allocation buffered writes highlighted in the scanner I/O
// benchmarks instead of allocating a new encoder per append.
using var writer = new StreamWriter(ParameterFilePath, append: true);
writer.WriteLine($"{length} {modulus.High} {modulus.Low} {primitiveRoot.High} {primitiveRoot.Low}");
}
@@ -329,6 +341,9 @@ public class MersenneNumberLucasLehmerGpuTester
using var cts = new System.Threading.CancellationTokenSource();
object sync = new();
// TODO: Move this Parallel.For to the shared low-overhead work scheduler once the NTT parameter
// generator integrates with the GPU-first pipeline so parameter scans reuse the same batching
// strategy measured fastest in the divisor-cycle benchmarks.
System.Threading.Tasks.Parallel.For(0, processorCount, (worker, state) =>
{
ulong k = (ulong)worker + 1UL;
@@ -377,6 +392,8 @@ public class MersenneNumberLucasLehmerGpuTester
private static ulong FindPrimitiveRoot(ulong modulus, ulong order)
{
ulong phi = modulus - 1UL;
// TODO: Reuse the divisor-cycle cache to factor phi via the precomputed small-prime windows once
// the lookup tables land so primitive root searches stop iterating over slow trial divisions.
var factors = Factorize(phi);
for (ulong g = 2UL; g < modulus; g++)
{
PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs
+11 -7
View File
@@ -23,11 +23,12 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO
GpuKernelType.Pow2Mod => gpuLease.Pow2ModOrderKernel,
_ => gpuLease.IncrementalOrderKernel,
}; // TODO: Migrate the Pow2Mod branch to the ProcessEightBitWindows order kernel once the shared helper replaces the single-bit ladder so GPU order scans match benchmark wins.
// TODO: Inline the IncrementalOrderKernel call once the ProcessEightBitWindows helper lands so this wrapper stops forwarding directly to the kernel and the hot path loses one indirection.
var foundBuffer = accelerator.Allocate1D<int>(1);
exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
try
{
var foundBuffer = accelerator.Allocate1D<int>(1); // TODO: Replace this allocation with the pooled device buffer from GpuOrderKernelBenchmarks so repeated scans reuse the pinned staging memory instead of allocating per run.
exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
try
{
UInt128 remaining;
int currentSize;
int found = 0;
@@ -37,9 +38,12 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
foundBuffer.MemSetToZero();
// Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + 1UL;
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
UInt128 q0 = twoP * kStart + 1UL;
// TODO: Replace these modulo recomputations with the divisor-cycle stepping deltas captured in
// MersenneDivisorCycleLengthGpuBenchmarks so each batch advances via cached residues from the
// single snapshot block instead of recalculating mod 10/8/5/3 in every iteration.
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
// Ensure device has small cycles table for early rejections in order kernels
PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs
+91 -87
View File
@@ -1,106 +1,110 @@
using System.Buffers;
using System.Runtime.InteropServices;
using System.Threading;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
public class MersenneNumberResidueGpuTester(bool useGpuOrder)
{
private readonly bool _useGpuOrder = useGpuOrder;
// TODO: Integrate MersenneDivisorCycles.Shared to consult cycle lengths for small q (<= 4M) and fast-reject
// candidates without launching heavy kernels. Consider a device-visible constant buffer per-accelerator
// via GpuKernelPool to avoid host round-trips.
private readonly bool _useGpuOrder = useGpuOrder;
// TODO: Integrate MersenneDivisorCycles.Shared to consult cycle lengths for small q (<= 4M) and fast-reject
// candidates without launching heavy kernels. Consider a device-visible constant buffer per-accelerator
// via GpuKernelPool to avoid host round-trips.
// TODO: When cycles are missing for larger q values, compute only the single required cycle on the device
// requested by the caller, skip queuing additional block generation, and keep the snapshot cache untouched.
// GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1.
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
// Ensure device has small cycles and primes tables for in-kernel lookup
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
int batchSize = GpuConstants.ScanBatchSize;
UInt128 kStart = 1UL;
UInt128 limit = maxK + UInt128.One;
byte last = lastIsSeven ? (byte)1 : (byte)0;
var kernel = gpuLease.Pow2ModKernel; // TODO: Swap this to the ProcessEightBitWindows kernel once GpuUInt128.Pow2Minus1Mod adopts the shared windowed helper measured fastest in GpuPow2ModBenchmarks.
twoP.Mod10_8_5_3(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
step10 = step10.Mod10();
// GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1.
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
// Ensure device has small cycles and primes tables for in-kernel lookup
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
int batchSize = GpuConstants.ScanBatchSize;
UInt128 kStart = 1UL;
UInt128 limit = maxK + UInt128.One;
byte last = lastIsSeven ? (byte)1 : (byte)0;
var kernel = gpuLease.Pow2ModKernel; // TODO: Swap this to the ProcessEightBitWindows kernel once GpuUInt128.Pow2Minus1Mod adopts the shared windowed helper measured fastest in GpuPow2ModBenchmarks.
twoP.Mod10_8_5_3(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
step10 = step10.Mod10();
GpuUInt128 twoPGpu = (GpuUInt128)twoP;
GpuUInt128 twoPGpu = (GpuUInt128)twoP;
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
// Device memory returned by Allocate1D is not guaranteed to be zeroed.
// Ensure the buffer starts with a known state so that lanes skipped by the
// kernel (for example due to early rejections) do not leave stale garbage
// values that would be interpreted as composite witnesses when copied back
// to the host. This mirrors the explicit zeroing performed by other GPU
// testers such as the order kernels.
orderBuffer.MemSetToZero();
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
UInt128 batchSize128 = (UInt128)batchSize;
UInt128 fullBatchStep = twoP * batchSize128;
UInt128 q = twoP * kStart + UInt128.One;
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
// Device memory returned by Allocate1D is not guaranteed to be zeroed.
// Ensure the buffer starts with a known state so that lanes skipped by the
// kernel (for example due to early rejections) do not leave stale garbage
// values that would be interpreted as composite witnesses when copied back
// to the host. This mirrors the explicit zeroing performed by other GPU
// testers such as the order kernels.
orderBuffer.MemSetToZero();
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
UInt128 batchSize128 = (UInt128)batchSize;
UInt128 fullBatchStep = twoP * batchSize128;
UInt128 q = twoP * kStart + UInt128.One;
// TODO: Replace the direct UInt128 multiply/add updates below with the residue-cycle stepping helper so that q advances reuse cached cycle increments rather than recomputing twoP multiples on every iteration.
try
{
while (kStart < limit && Volatile.Read(ref isPrime))
{
UInt128 remaining = limit - kStart;
int currentSize = remaining > batchSize128 ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize);
ref ulong ordersRef = ref MemoryMarshal.GetReference(orders);
try
{
while (kStart < limit && Volatile.Read(ref isPrime))
{
UInt128 remaining = limit - kStart;
int currentSize = remaining > batchSize128 ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize);
ref ulong ordersRef = ref MemoryMarshal.GetReference(orders);
q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL,
kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2);
kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL,
kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2);
accelerator.Synchronize();
orderBuffer.View.CopyToCPU(ref ordersRef, currentSize);
if (!Volatile.Read(ref isPrime))
{
break;
}
accelerator.Synchronize();
orderBuffer.View.CopyToCPU(ref ordersRef, currentSize);
if (!Volatile.Read(ref isPrime))
{
break;
}
bool compositeFound = false;
for (int i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
compositeFound = true;
break;
}
}
bool compositeFound = false;
for (int i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
compositeFound = true;
break;
}
}
if (compositeFound)
{
break;
}
if (compositeFound)
{
break;
}
UInt128 processed = (UInt128)currentSize;
kStart += processed;
if (kStart < limit)
{
if (processed == batchSize128)
{
q += fullBatchStep;
}
else
{
q += twoP * processed;
}
}
}
}
finally
{
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
gpuLease.Dispose();
}
}
UInt128 processed = (UInt128)currentSize;
kStart += processed;
if (kStart < limit)
{
if (processed == batchSize128)
{
q += fullBatchStep;
}
else
{
q += twoP * processed; // TODO: Swap this multiplication for the shared cycle stepping helper once residue cycles are mandatory so we can reuse the cached remainder ladder instead of performing UInt128 multiply operations.
}
}
}
}
finally
{
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
gpuLease.Dispose();
}
}
}
PerfectNumbers.Core/Gpu/NttGpuMath.cs
+21 -14
View File
@@ -29,23 +29,23 @@ public static class NttGpuMath
public static NttBackend GpuTransformBackend { get; set; } = NttBackend.Reference;
// Controls modular reduction strategy inside staged NTT paths.
public static ModReductionMode ReductionMode { get; set; } = ModReductionMode.Auto;
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<GpuUInt128>, GpuUInt128>> MulKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int, int, int, ArrayView<GpuUInt128>, GpuUInt128>> StageKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, GpuUInt128, GpuUInt128>> ScaleKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int, int, int, ArrayView<GpuUInt128>, ulong, ulong>> StageMontKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int, int, int, ArrayView<GpuUInt128>, ulong, ulong, ulong, ulong>> StageBarrett128KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, GpuUInt128, ulong, ulong, ulong, ulong>> ScaleBarrett128KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong, ulong, ulong>> SquareBarrett128KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong, ulong>> ToMont64KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong>> FromMont64KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong>> SquareMont64KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong, ulong>> ScaleMont64KernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<GpuUInt128>, GpuUInt128>> MulKernelCache = new(); // TODO: Replace this concurrent cache with the prewarmed accelerator-indexed tables from GpuModularArithmeticBenchmarks so kernel launches avoid dictionary lookups once the kernels are baked during startup.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int, int, int, ArrayView<GpuUInt128>, GpuUInt128>> StageKernelCache = new(); // TODO: Same as above – materialize staged kernels during startup so we can drop ConcurrentDictionary usage entirely per the benchmark guidance.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, GpuUInt128, GpuUInt128>> ScaleKernelCache = new(); // TODO: Promote to the static kernel table initialized alongside the GpuModularArithmeticBenchmarks fast path.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int, int, int, ArrayView<GpuUInt128>, ulong, ulong>> StageMontKernelCache = new(); // TODO: Inline the Montgomery kernels into the startup table instead of using a concurrent cache now that we no longer mutate state at runtime.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int, int, int, ArrayView<GpuUInt128>, ulong, ulong, ulong, ulong>> StageBarrett128KernelCache = new(); // TODO: Same plan – reuse the precomputed kernels measured fastest in MontgomeryMultiplyBenchmarks rather than looking them up dynamically.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, GpuUInt128, ulong, ulong, ulong, ulong>> ScaleBarrett128KernelCache = new(); // TODO: Collapse into the startup kernel array once Barrett128 constants are preloaded.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong, ulong, ulong>> SquareBarrett128KernelCache = new(); // TODO: Fold these kernels into the same startup table to avoid concurrent access overhead highlighted in GpuModularArithmeticBenchmarks.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong, ulong>> ToMont64KernelCache = new(); // TODO: Prebind and reuse the Montgomery conversions from MontgomeryMultiplyBenchmarks instead of storing them in a concurrent cache.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong>> FromMont64KernelCache = new(); // TODO: As above – drop ConcurrentDictionary once startup prewarms every kernel.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong>> SquareMont64KernelCache = new(); // TODO: Inline into the static kernel table established during initialization.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ulong, ulong, ulong>> ScaleMont64KernelCache = new(); // TODO: Move to the startup kernel table so Montgomery-scaled launches skip concurrent lookups.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<GpuUInt128>, int, GpuUInt128, GpuUInt128>> ForwardKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<GpuUInt128>, int, GpuUInt128, GpuUInt128>> ForwardKernelCache = new(); // TODO: Warm these forward kernels during initialization and store them in a plain array so the LucasLehmerGpuBenchmarks launch costs disappear.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<GpuUInt128>, int, GpuUInt128, GpuUInt128, GpuUInt128>> InverseKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<GpuUInt128>, int, GpuUInt128, GpuUInt128, GpuUInt128>> InverseKernelCache = new(); // TODO: Same as above – replace with startup-prepared tables per LucasLehmerGpuBenchmarks guidance.
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int>> BitReverseKernelCache = new();
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int>> BitReverseKernelCache = new(); // TODO: Precompute bit-reversal kernels during boot so we can remove the concurrent lookup overhead flagged in the benchmarks.
private sealed class SquareCacheEntry : IDisposable
{
@@ -111,6 +111,9 @@ public static class NttGpuMath
int half = len >> 1;
StageOffsets[stage - 1] = offset;
UInt128 exp = (modValue - 1UL) / (ulong)len;
// TODO: Pull these stage roots from the precomputed tables measured fastest in
// Pow2MontgomeryModCycleComputationBenchmarks instead of recomputing ModPow for
// every stage.
var wLen = new GpuUInt128(Root);
wLen.ModPow((ulong)exp, modulus);
var w = new GpuUInt128(0UL, 1UL);
@@ -119,6 +122,8 @@ public static class NttGpuMath
forward[offset + j] = w;
w.MulMod(wLen, modulus);
}
// TODO: Reuse the same precomputed tables for inverse roots so we do not repeat the
// expensive ModPow calls highlighted as bottlenecks in Pow2MontgomeryModCycleComputationBenchmarks.
var wLenInv = new GpuUInt128(RootInv);
wLenInv.ModPow((ulong)exp, modulus);
var wi = new GpuUInt128(0UL, 1UL);
@@ -140,6 +145,8 @@ public static class NttGpuMath
{
UseMontgomery64 = true;
ModulusLow = modulus.Low;
// TODO: Preload MontNPrime64/MontR values using the MontgomeryMultiplyBenchmarks helper so
// startup reuses the fastest constant generation routine instead of recomputing here.
MontNPrime64 = ComputeMontgomeryNPrime64(ModulusLow);
var mBI = new BigInteger(ModulusLow);
MontRMod64 = (ulong)((BigInteger.One << 64) % mBI);