From fddec17b3f9575d186369cded5ed35dd34b54243 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Fri, 17 Oct 2025 06:53:54 +0200 Subject: [PATCH] Various try-out changes --- PerfectNumbers.Core/Gpu/GpuKernelLease.cs | 707 ++++++++++++++++++ PerfectNumbers.Core/Gpu/GpuKernelPool.cs | 702 ----------------- .../PrimeOrderCalculator.Cpu.cs | 29 +- .../PrimeOrderCalculator.Gpu.cs | 156 ++-- .../PrimeOrderCalculator.Shared.cs | 38 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 2 +- 6 files changed, 813 insertions(+), 821 deletions(-) create mode 100644 PerfectNumbers.Core/Gpu/GpuKernelLease.cs diff --git a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs new file mode 100644 index 0000000..2318f6d --- /dev/null +++ b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs @@ -0,0 +1,707 @@ +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Runtime; +using static PerfectNumbers.Core.Gpu.GpuContextPool; + +namespace PerfectNumbers.Core.Gpu; + +public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels) +{ + private bool disposedValue; + private IDisposable? _limiter = limiter; + private GpuContextLease? _gpu = gpu; + private KernelContainer _kernels = kernels; + private AcceleratorStream? _stream; + private object? _executionLock = gpu.ExecutionLock; + + public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released"); + + public AcceleratorStream Stream + { + get + { + if (_stream is { } stream) + { + return stream; + } + + // It feels that we can combine these into 1 return statement, but that increases the max stack frame size to 3 from 2 + stream = Accelerator.CreateStream(); + _stream = stream; + return stream; + } + } + + public readonly ExecutionScope EnterExecutionScope() => new ExecutionScope(_executionLock); + + public readonly Action, ArrayView> OrderKernel + { + get + { + var accel = Accelerator; // avoid capturing 'this' in lambda + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Order, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView>>(); + }); + } + } + } + + public readonly Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> Pow2ModKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + }); + } + } + } + + public readonly Action, ArrayView1D> IncrementalKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Incremental, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D>>(); + }); + } + } + } + + public readonly Action, ArrayView1D> IncrementalOrderKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalOrderKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D>>(); + }); + } + } + } + + public readonly Action, ArrayView1D> Pow2ModOrderKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModOrderKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D>>(); + }); + } + } + } + + public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.SmallPrimeFactor, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + }); + } + } + } + + private static void SmallPrimeFactorKernelScan( + Index1D index, + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D primeSlots, + ArrayView1D exponentSlots, + ArrayView1D countSlot, + ArrayView1D remainingSlot) + { + if (index != 0) + { + return; + } + + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + ulong remaining = value; + int count = 0; + int capacity = (int)primeSlots.Length; + + for (int i = 0; i < primeCount && remaining > 1UL && count < capacity; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + primeSlots[count] = primeValue; + exponentSlots[count] = exponent; + count++; + } + + countSlot[0] = count; + remainingSlot[0] = remaining; + } + + // TODO: Plumb the small-cycles device buffer into all kernels that can benefit + // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. + + private static void IncrementalKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, + ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView orders, + ArrayView1D smallCycles) + { + ulong idx = (ulong)index.X; + // TODO: Replace these `%` computations with the precomputed Mod3/Mod5 tables so GPU kernels reuse cached residues + // instead of performing modulo operations that the benchmarks showed slower on wide sweeps. + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + // residue automaton + ulong step10 = (exponent.Mod10() << 1).Mod10(); + ulong r10 = q0m10 + (step10 * idx).Mod10(); + r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5 + if (shouldCheck) + { + ulong step8 = ((exponent & 7UL) << 1) & 7UL; + ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked + // bitmask-based implementation for GPU workloads. + ulong step3 = ((exponent % 3UL) << 1) % 3UL; + ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong step5Loc = ((exponent % 5UL) << 1) % 5UL; + ulong r5 = q0m5 + (step5Loc * idxMod5); if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + + if (!shouldCheck) + { + orders[index] = 0UL; + return; + } + + // TODO: Is this expected that the value of kStart is modified? + GpuUInt128 k = kStart + (GpuUInt128)idx; + GpuUInt128 q = twoP; + q.Mul64(k); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + + // Small-cycles in-kernel early rejection from device table + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + orders[index] = 0UL; + return; + } + } + + // TODO: q will be lowered by GpuUInt128.One here. Is this expected? Maybe that's why we have issues with the incremental kernel? Review other places where we use GpuUInt128 operators. + GpuUInt128 phi = q - GpuUInt128.One; + // TODO: Do we have the check here because we don't support high values or this is expected result? + if (phi.High != 0UL) + { + orders[index] = 0UL; + return; + } + + ulong phi64 = phi.Low; + // TODO: Replace these Pow2Mod calls with the ProcessEightBitWindows helper when the shared windowed + // scalar implementation lands; benchmarks showed the windowed kernel trimming per-divisor runtime by ~2.4×. + if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + ulong div = FastDiv64Gpu(phi64, exponent, divMul); + GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + orders[index] = 1UL; + } + + private static void OrderKernelScan(Index1D index, ulong exponent, ulong divMul, ArrayView qs, ArrayView orders) + { + GpuUInt128 q = qs[index]; + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + GpuUInt128 phi = q - GpuUInt128.One; + if (phi.High != 0UL) + { + orders[index] = 0UL; + return; + } + + ulong phi64 = phi.Low; + // TODO: Once the ProcessEightBitWindows helper is available, switch this order kernel to that faster + // Pow2Mod variant so cycle checks inherit the same gains observed in GpuPow2ModBenchmarks. + GpuUInt128 pow = GpuUInt128.Pow2Mod(phi64, in readOnlyQ); + if (pow != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + ulong div = FastDiv64Gpu(phi64, exponent, divMul); + GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + orders[index] = exponent; + } + + private void Dispose(bool disposing) + { + if (!disposedValue) + { + if (disposing) + { + _stream?.Dispose(); + _stream = null; + + _gpu?.Dispose(); + _gpu = null; + + _executionLock = null; + + _limiter?.Dispose(); + _limiter = null; + } + + // free unmanaged resources (unmanaged objects) and override finalizer + // set large fields to null + disposedValue = true; + } + } + + private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, + ResidueAutomatonArgs ra, ArrayView orders, + ArrayView1D smallCycles, + ArrayView1D smallPrimesLastOne, + ArrayView1D smallPrimesLastSeven, + ArrayView1D smallPrimesPow2LastOne, + ArrayView1D smallPrimesPow2LastSeven) + { + ulong idx = (ulong)index.X; + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + orders[index] = 0UL; + return; + } + + // TODO: kStart is modified after this. Is this expected? + kStart.Add(idx); + GpuUInt128 q = twoP; + q.Mul64(kStart); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + // cycle should always be initialized if we're within array limit in production code + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + orders[index] = 0UL; + return; + } + } + // TODO: Swap Pow2Minus1Mod for the eight-bit window helper once the scalar version switches; + // benchmarks show the windowed variant cuts large-divisor scans from ~51 µs to ~21 µs. + if (GpuUInt128.Pow2Minus1Mod(exponent, in readOnlyQ) != GpuUInt128.Zero) + { + orders[index] = 0UL; + return; + } + + ArrayView1D primes = lastIsSeven != 0 ? smallPrimesLastSeven : smallPrimesLastOne; + ArrayView1D primesPow2 = lastIsSeven != 0 ? smallPrimesPow2LastSeven : smallPrimesPow2LastOne; + int primesLen = (int)primes.Length; + for (int i = 0; i < primesLen; i++) + { + ulong square = primesPow2[i]; + if (new GpuUInt128(0UL, square) > q) + { + break; + } + ulong prime = primes[i]; + if (Mod128By64(q, prime) == 0UL) + { + orders[index] = 0UL; + return; + } + } + + orders[index] = exponent; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong FastDiv64Gpu(ulong value, ulong divisor, ulong mul) + { + ulong quotient = GpuUInt128.MulHigh(value, mul); + GpuUInt128 remainder = new(0UL, value); + GpuUInt128 product = new(GpuUInt128.MulHigh(quotient, divisor), quotient * divisor); + remainder.Sub(product); + + if (remainder.High != 0UL || remainder.Low >= divisor) + { + quotient++; + } + + return quotient; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Mod128By64(GpuUInt128 value, ulong modulus) + { + ulong rem = 0UL; + ulong part = value.High; + for (int i = 0; i < 64; i++) + { + rem = (rem << 1) | (part >> 63); + if (rem >= modulus) + { + rem -= modulus; + } + part <<= 1; + } + + part = value.Low; + for (int i = 0; i < 64; i++) + { + rem = (rem << 1) | (part >> 63); + if (rem >= modulus) + { + rem -= modulus; + } + part <<= 1; + } + + return rem; + } + + // Device-friendly small cycle length for 2 mod divisor + private static ulong CalculateCycleLengthSmall(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + return 1UL; + + ulong order = 1UL, pow = 2UL; + while (pow != 1UL) + { + pow <<= 1; + if (pow > divisor) + pow -= divisor; + + order++; + } + + return order; + } + + private static void IncrementalOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, + ResidueAutomatonArgs ra, ArrayView found, ArrayView1D smallCycles) + { + ulong idx = (ulong)index.X; + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + return; + } + + GpuUInt128 k = kStart + (GpuUInt128)idx; + GpuUInt128 q = twoP; + q.Mul64(k); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + + // Small-cycles in-kernel early rejection from device table + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + return; + } + } + GpuUInt128 phi = q - GpuUInt128.One; + if (phi.High != 0UL) + { + return; + } + + ulong phi64 = phi.Low; + // TODO: Upgrade this pow2mod order kernel to the ProcessEightBitWindows helper once available so GPU residue + // scans avoid the single-bit ladder that benchmarks found to be 2.3× slower on large exponents. + if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) + { + return; + } + + GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) + { + return; + } + + ulong div = FastDiv64Gpu(phi64, exponent, divMul); + GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) + { + return; + } + + Atomic.Or(ref found[0], 1); + } + + private static void Pow2ModOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, + ResidueAutomatonArgs ra, ArrayView found, ArrayView1D smallCycles) + { + ulong idx = (ulong)index.X; + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + return; + } + + GpuUInt128 k = kStart + (GpuUInt128)idx; + GpuUInt128 q = twoP; + q.Mul64(k); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + // Small-cycles in-kernel early rejection from device table + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + return; + } + } + // TODO: Replace this Pow2Mod check with the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; + // residue order scans will then benefit from the same 2× speedup the benchmarked windowed kernel delivered. + if (GpuUInt128.Pow2Mod(exponent, in readOnlyQ) != GpuUInt128.One) + { + return; + } + + Atomic.Or(ref found[0], 1); + } + + public void Dispose() + { + // Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method + Dispose(disposing: true); + // Enabled this when the struct is converted to a class + // GC.SuppressFinalize(this); + } + + public readonly struct ExecutionScope + { + private readonly object? _lock; + + public ExecutionScope(object? sync) + { + _lock = sync; + if (_lock is not null) + { + // Monitor.Enter(_lock); + } + } + + public void Dispose() + { + if (_lock is not null) + { + // Monitor.Exit(_lock); + } + } + } +} diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index cf03ebb..2e63cd0 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -1,5 +1,4 @@ using System.Collections.Concurrent; -using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; using static PerfectNumbers.Core.Gpu.GpuContextPool; @@ -130,707 +129,6 @@ public readonly struct SmallPrimeFactorScratch } } -public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels) -{ - private bool disposedValue; - private IDisposable? _limiter = limiter; - private GpuContextLease? _gpu = gpu; - private KernelContainer _kernels = kernels; - private AcceleratorStream? _stream; - private object? _executionLock = gpu.ExecutionLock; - - public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released"); - - public AcceleratorStream Stream - { - get - { - if (_stream is { } stream) - { - return stream; - } - - // It feels that we can combine these into 1 return statement, but that increases the max stack frame size to 3 from 2 - stream = Accelerator.CreateStream(); - _stream = stream; - return stream; - } - } - - public readonly ExecutionScope EnterExecutionScope() => new ExecutionScope(_executionLock); - - public readonly Action, ArrayView> OrderKernel - { - get - { - var accel = Accelerator; // avoid capturing 'this' in lambda - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Order, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView>>(); - }); - } - } - } - - public readonly Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> Pow2ModKernel - { - get - { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); - }); - } - } - } - - public readonly Action, ArrayView1D> IncrementalKernel - { - get - { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Incremental, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D>>(); - }); - } - } - } - - public readonly Action, ArrayView1D> IncrementalOrderKernel - { - get - { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalOrderKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D>>(); - }); - } - } - } - - public readonly Action, ArrayView1D> Pow2ModOrderKernel - { - get - { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModOrderKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D>>(); - }); - } - } - } - - public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernel - { - get - { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.SmallPrimeFactor, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); - }); - } - } - } - - private static void SmallPrimeFactorKernelScan( - Index1D index, - ulong value, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D primeSlots, - ArrayView1D exponentSlots, - ArrayView1D countSlot, - ArrayView1D remainingSlot) - { - if (index != 0) - { - return; - } - - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - ulong remaining = value; - int count = 0; - int capacity = (int)primeSlots.Length; - - for (int i = 0; i < primeCount && remaining > 1UL && count < capacity; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate == 0U || primeCandidate > effectiveLimit) - { - break; - } - - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remaining) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == 0UL); - - primeSlots[count] = primeValue; - exponentSlots[count] = exponent; - count++; - } - - countSlot[0] = count; - remainingSlot[0] = remaining; - } - - // TODO: Plumb the small-cycles device buffer into all kernels that can benefit - // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. - - private static void IncrementalKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, - ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView orders, - ArrayView1D smallCycles) - { - ulong idx = (ulong)index.X; - // TODO: Replace these `%` computations with the precomputed Mod3/Mod5 tables so GPU kernels reuse cached residues - // instead of performing modulo operations that the benchmarks showed slower on wide sweeps. - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - // residue automaton - ulong step10 = (exponent.Mod10() << 1).Mod10(); - ulong r10 = q0m10 + (step10 * idx).Mod10(); - r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5 - if (shouldCheck) - { - ulong step8 = ((exponent & 7UL) << 1) & 7UL; - ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked - // bitmask-based implementation for GPU workloads. - ulong step3 = ((exponent % 3UL) << 1) % 3UL; - ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong step5Loc = ((exponent % 5UL) << 1) % 5UL; - ulong r5 = q0m5 + (step5Loc * idxMod5); if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - - if (!shouldCheck) - { - orders[index] = 0UL; - return; - } - - // TODO: Is this expected that the value of kStart is modified? - GpuUInt128 k = kStart + (GpuUInt128)idx; - GpuUInt128 q = twoP; - q.Mul64(k); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - - // Small-cycles in-kernel early rejection from device table - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - orders[index] = 0UL; - return; - } - } - - // TODO: q will be lowered by GpuUInt128.One here. Is this expected? Maybe that's why we have issues with the incremental kernel? Review other places where we use GpuUInt128 operators. - GpuUInt128 phi = q - GpuUInt128.One; - // TODO: Do we have the check here because we don't support high values or this is expected result? - if (phi.High != 0UL) - { - orders[index] = 0UL; - return; - } - - ulong phi64 = phi.Low; - // TODO: Replace these Pow2Mod calls with the ProcessEightBitWindows helper when the shared windowed - // scalar implementation lands; benchmarks showed the windowed kernel trimming per-divisor runtime by ~2.4×. - if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - ulong div = FastDiv64Gpu(phi64, exponent, divMul); - GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - orders[index] = 1UL; - } - - private static void OrderKernelScan(Index1D index, ulong exponent, ulong divMul, ArrayView qs, ArrayView orders) - { - GpuUInt128 q = qs[index]; - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - GpuUInt128 phi = q - GpuUInt128.One; - if (phi.High != 0UL) - { - orders[index] = 0UL; - return; - } - - ulong phi64 = phi.Low; - // TODO: Once the ProcessEightBitWindows helper is available, switch this order kernel to that faster - // Pow2Mod variant so cycle checks inherit the same gains observed in GpuPow2ModBenchmarks. - GpuUInt128 pow = GpuUInt128.Pow2Mod(phi64, in readOnlyQ); - if (pow != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - ulong div = FastDiv64Gpu(phi64, exponent, divMul); - GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - orders[index] = exponent; - } - - private void Dispose(bool disposing) - { - if (!disposedValue) - { - if (disposing) - { - _stream?.Dispose(); - _stream = null; - - _gpu?.Dispose(); - _gpu = null; - - _executionLock = null; - - _limiter?.Dispose(); - _limiter = null; - } - - // free unmanaged resources (unmanaged objects) and override finalizer - // set large fields to null - disposedValue = true; - } - } - - private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, - ResidueAutomatonArgs ra, ArrayView orders, - ArrayView1D smallCycles, - ArrayView1D smallPrimesLastOne, - ArrayView1D smallPrimesLastSeven, - ArrayView1D smallPrimesPow2LastOne, - ArrayView1D smallPrimesPow2LastSeven) - { - ulong idx = (ulong)index.X; - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - orders[index] = 0UL; - return; - } - - // TODO: kStart is modified after this. Is this expected? - kStart.Add(idx); - GpuUInt128 q = twoP; - q.Mul64(kStart); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - // cycle should always be initialized if we're within array limit in production code - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - orders[index] = 0UL; - return; - } - } - // TODO: Swap Pow2Minus1Mod for the eight-bit window helper once the scalar version switches; - // benchmarks show the windowed variant cuts large-divisor scans from ~51 µs to ~21 µs. - if (GpuUInt128.Pow2Minus1Mod(exponent, in readOnlyQ) != GpuUInt128.Zero) - { - orders[index] = 0UL; - return; - } - - ArrayView1D primes = lastIsSeven != 0 ? smallPrimesLastSeven : smallPrimesLastOne; - ArrayView1D primesPow2 = lastIsSeven != 0 ? smallPrimesPow2LastSeven : smallPrimesPow2LastOne; - int primesLen = (int)primes.Length; - for (int i = 0; i < primesLen; i++) - { - ulong square = primesPow2[i]; - if (new GpuUInt128(0UL, square) > q) - { - break; - } - ulong prime = primes[i]; - if (Mod128By64(q, prime) == 0UL) - { - orders[index] = 0UL; - return; - } - } - - orders[index] = exponent; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong FastDiv64Gpu(ulong value, ulong divisor, ulong mul) - { - ulong quotient = GpuUInt128.MulHigh(value, mul); - GpuUInt128 remainder = new(0UL, value); - GpuUInt128 product = new(GpuUInt128.MulHigh(quotient, divisor), quotient * divisor); - remainder.Sub(product); - - if (remainder.High != 0UL || remainder.Low >= divisor) - { - quotient++; - } - - return quotient; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong Mod128By64(GpuUInt128 value, ulong modulus) - { - ulong rem = 0UL; - ulong part = value.High; - for (int i = 0; i < 64; i++) - { - rem = (rem << 1) | (part >> 63); - if (rem >= modulus) - { - rem -= modulus; - } - part <<= 1; - } - - part = value.Low; - for (int i = 0; i < 64; i++) - { - rem = (rem << 1) | (part >> 63); - if (rem >= modulus) - { - rem -= modulus; - } - part <<= 1; - } - - return rem; - } - - // Device-friendly small cycle length for 2 mod divisor - private static ulong CalculateCycleLengthSmall(ulong divisor) - { - if ((divisor & (divisor - 1UL)) == 0UL) - return 1UL; - - ulong order = 1UL, pow = 2UL; - while (pow != 1UL) - { - pow <<= 1; - if (pow > divisor) - pow -= divisor; - - order++; - } - - return order; - } - - private static void IncrementalOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, - ResidueAutomatonArgs ra, ArrayView found, ArrayView1D smallCycles) - { - ulong idx = (ulong)index.X; - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - return; - } - - GpuUInt128 k = kStart + (GpuUInt128)idx; - GpuUInt128 q = twoP; - q.Mul64(k); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - - // Small-cycles in-kernel early rejection from device table - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - return; - } - } - GpuUInt128 phi = q - GpuUInt128.One; - if (phi.High != 0UL) - { - return; - } - - ulong phi64 = phi.Low; - // TODO: Upgrade this pow2mod order kernel to the ProcessEightBitWindows helper once available so GPU residue - // scans avoid the single-bit ladder that benchmarks found to be 2.3× slower on large exponents. - if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) - { - return; - } - - GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) - { - return; - } - - ulong div = FastDiv64Gpu(phi64, exponent, divMul); - GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) - { - return; - } - - Atomic.Or(ref found[0], 1); - } - - private static void Pow2ModOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, - ResidueAutomatonArgs ra, ArrayView found, ArrayView1D smallCycles) - { - ulong idx = (ulong)index.X; - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - return; - } - - GpuUInt128 k = kStart + (GpuUInt128)idx; - GpuUInt128 q = twoP; - q.Mul64(k); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - // Small-cycles in-kernel early rejection from device table - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - return; - } - } - // TODO: Replace this Pow2Mod check with the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; - // residue order scans will then benefit from the same 2× speedup the benchmarked windowed kernel delivered. - if (GpuUInt128.Pow2Mod(exponent, in readOnlyQ) != GpuUInt128.One) - { - return; - } - - Atomic.Or(ref found[0], 1); - } - - public void Dispose() - { - // Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method - Dispose(disposing: true); - // Enabled this when the struct is converted to a class - // GC.SuppressFinalize(this); - } - - public readonly struct ExecutionScope - { - private readonly object? _lock; - - public ExecutionScope(object? sync) - { - _lock = sync; - if (_lock is not null) - { - Monitor.Enter(_lock); - } - } - - public void Dispose() - { - if (_lock is not null) - { - Monitor.Exit(_lock); - } - } - } -} - public class GpuKernelPool { private static readonly ConcurrentDictionary 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. diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 50357bd..41c7626 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -733,16 +733,17 @@ internal static partial class PrimeOrderCalculator counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(FactorSlotCount, 8)); counts.Clear(); - bool gpuPopulated = TryPopulateSmallPrimeFactorsGpu(value, limit, counts, out remaining); - if (!gpuPopulated) - { - counts.Clear(); + // bool gpuPopulated = TryPopulateSmallPrimeFactorsGpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining); + // if (!gpuPopulated) + // { + // throw new InvalidOperationException($"GPU didn't populate factors for {value}"); + // counts.Clear(); remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); - } + // } - int dictionaryCount = counts.Count; - if (dictionaryCount <= FactorSlotCount) - { + // int dictionaryCount = counts.Count; + // if (dictionaryCount <= FactorSlotCount) + // { int copyIndex = 0; foreach (KeyValuePair entry in counts) { @@ -754,12 +755,12 @@ internal static partial class PrimeOrderCalculator factorCount = copyIndex; ThreadStaticPools.ReturnUlongIntDictionary(counts); counts = null; - } - else - { - useDictionary = true; - factorCount = dictionaryCount; - } + // } + // else + // { + // useDictionary = true; + // factorCount = dictionaryCount; + // } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index f21d726..1073e91 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -9,100 +9,86 @@ namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; + private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; - private const int GpuSmallPrimeFactorSlots = 64; + private const int GpuSmallPrimeFactorSlots = 64; - private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) - { - ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); - int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); - Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - primeBuffer.Clear(); - exponentBuffer.Clear(); - remaining = value; + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Span primeBuffer, Span exponentBuffer, out int factorCount, out ulong remaining) + { + // primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); + // exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); + // Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + // Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + // primeBuffer.Clear(); + // exponentBuffer.Clear(); + remaining = value; - try - { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - try - { - var execution = lease.EnterExecutionScope(); - try - { - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; + // try + // { + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + try + { + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; - SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator); - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); - scratch.Clear(); + SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator); + SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); + // scratch.Clear(); - var kernel = lease.SmallPrimeFactorKernel; - kernel( - stream, - 1, - value, - limit, - tables.PrimesView, - tables.SquaresView, - tables.Count, - scratch.PrimeSlots.View, - scratch.ExponentSlots.View, - scratch.CountSlot.View, - scratch.RemainingSlot.View); + var kernel = lease.SmallPrimeFactorKernel; + kernel( + stream, + 1, + value, + limit, + tables.PrimesView, + tables.SquaresView, + tables.Count, + scratch.PrimeSlots.View, + scratch.ExponentSlots.View, + scratch.CountSlot.View, + scratch.RemainingSlot.View); - stream.Synchronize(); + stream.Synchronize(); - int factorCount = 0; - scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); - factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); + factorCount = 0; + scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); + factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); - if (factorCount > 0) - { - scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount); - scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount); - } - scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); + if (factorCount > 0) + { + scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount); + scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount); + } - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - if (primeValue == 0UL || exponent == 0) - { - continue; - } + scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); - counts[primeValue] = exponent; - } + // for (int i = 0; i < factorCount; i++) + // { + // ulong primeValue = primeBuffer[i]; + // int exponent = exponentBuffer[i]; + // // This will never happen in production code + // // if (primeValue == 0UL || exponent == 0) + // // { + // // continue; + // // } - return true; - } - finally - { - execution.Dispose(); - } - } - finally - { - lease.Dispose(); - } - } - catch (CLException) - { - remaining = value; - return false; - } - catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) - { - remaining = value; - return false; - } - finally - { - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - } - } + // counts[primeValue] = exponent; + // } + + return true; + } + finally + { + execution.Dispose(); + lease.Dispose(); + } + // } + // finally + // { + // ThreadStaticPools.UlongPool.Return(primeBufferArray); + // ThreadStaticPools.IntPool.Return(exponentBufferArray); + // } + } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs index d586305..d385fb8 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -606,27 +606,27 @@ internal static partial class PrimeOrderCalculator uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; UInt128 remaining; - if (IsGpuPow2Allowed && value <= ulong.MaxValue) - { - Dictionary narrowCounts = new(capacity: 8); - if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining)) - { - foreach (KeyValuePair entry in narrowCounts) - { - counts[(UInt128)entry.Key] = entry.Value; - } + // if (IsGpuPow2Allowed && value <= ulong.MaxValue) + // { + // Dictionary narrowCounts = new(capacity: 8); + // if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining)) + // { + // foreach (KeyValuePair entry in narrowCounts) + // { + // counts[(UInt128)entry.Key] = entry.Value; + // } - remaining = (UInt128)narrowRemaining; - } - else - { - remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - } - } - else - { + // remaining = (UInt128)narrowRemaining; + // } + // else + // { + // remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); + // } + // } + // else + // { remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - } + // } List pending = new(); if (remaining > UInt128.One) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index d6ed485..bd8a473 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -10,7 +10,7 @@ internal static partial class PrimeOrderCalculator internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) { - public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 10_240, maxPowChecks: 24, PrimeOrderMode.Heuristic); + public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24, maxPowChecks: 24, PrimeOrderMode.Heuristic); public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); public readonly uint SmallFactorLimit = smallFactorLimit;