diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 48e8e27..5240b52 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -64,6 +64,7 @@ internal static class Program GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); GpuContextPool.WarmUpPool(threadCount); PrimeTester.WarmUpGpuKernels(gpuPrimeThreads); + GpuScratchBufferPool.WarmUp(gpuPrimeThreads << 4, 64, 1024); Console.WriteLine("Starting up threads..."); _ = UnboundedTaskScheduler.Instance; int blockSize = Math.Max(1, _cliArguments.BlockSize); diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index 92cb06b..4018959 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -330,26 +330,26 @@ public class GpuKernelCompilationTests private static void CompileMersenneDivisorCyclesKernel(Accelerator accelerator) { - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView1D>(MersenneDivisorCyclesKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel(MersenneDivisorCyclesKernel); _ = kernel; } private static void CompilePrimeTesterSmallPrimeKernel(Accelerator accelerator) { - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterSmallPrimeKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel(PrimeTesterSmallPrimeKernel); _ = kernel; } private static void CompilePrimeTesterSharesFactorKernel(Accelerator accelerator) { - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterSharesFactorKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel(PrimeTesterSharesFactorKernel); _ = kernel; } private static void CompileGpuKernelPoolKernels(Accelerator accelerator) { _ = accelerator; - var lease = GpuKernelPool.GetKernel(useGpuOrder: false); + var lease = GpuKernelPool.GetKernel(); try { _ = lease.OrderKernel; @@ -363,7 +363,7 @@ public class GpuKernelCompilationTests lease.Dispose(); } - var orderLease = GpuKernelPool.GetKernel(useGpuOrder: true); + var orderLease = GpuKernelPool.GetKernel(); try { _ = orderLease.OrderKernel; diff --git a/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs b/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs index fa078b0..a250fd5 100644 --- a/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs @@ -150,7 +150,7 @@ public class NttGpuMathTests var modulus = new GpuUInt128(0UL, 17UL); var primitiveRoot = new GpuUInt128(0UL, 3UL); - GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent(); + GpuContextLease gpu = GpuContextPool.Rent(); try { var accelerator = gpu.Accelerator; @@ -345,7 +345,7 @@ public class NttGpuMathTests values[1] = new GpuUInt128(0UL, 2UL); var modulus = new GpuUInt128(0UL, 17UL); var primitiveRoot = new GpuUInt128(0UL, 3UL); - GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent(); + GpuContextLease gpu = GpuContextPool.Rent(); try { var buffer = gpu.Accelerator.Allocate1D(values.Length); diff --git a/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs b/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs index 57fe816..50440ca 100644 --- a/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs @@ -33,7 +33,7 @@ public class NttMontgomeryTests NttGpuMath.Square(cpu, modulus, primitiveRoot); // Device execution - GpuContextPool.GpuContextLease context = GpuContextPool.Rent(); + GpuContextLease context = GpuContextPool.Rent(); try { var accelerator = context.Accelerator; diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 47521b4..c6f1309 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -226,9 +226,10 @@ public sealed class DivisorCycleCache private void ComputeCyclesGpuCore(ReadOnlySpan divisors, Span destination) { int length = divisors.Length; - var gpuLease = GpuKernelPool.GetKernel(useGpuOrder: true); + var gpuLease = GpuKernelPool.GetKernel(); - Accelerator accelerator = gpuLease.Accelerator; + Accelerator accelerator = gpuLease.Accelerator; + var stream = gpuLease.Stream; var kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel); using MemoryBuffer1D divisorBuffer = accelerator.Allocate1D(length); @@ -258,7 +259,7 @@ public sealed class DivisorCycleCache : new Span(rentedStatus = ArrayPool.Shared.Rent(length), 0, length); ref ulong divisorRef = ref MemoryMarshal.GetReference(divisors); - divisorBuffer.View.CopyFromCPU(ref divisorRef, length); + divisorBuffer.View.CopyFromCPU(stream, ref divisorRef, length); for (int i = 0; i < length; i++) { @@ -275,18 +276,18 @@ public sealed class DivisorCycleCache statusSpan[i] = ByteZero; } - powBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(powSpan), length); - orderBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(orderSpan), length); - resultBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(resultSpan), length); - statusBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(statusSpan), length); + powBuffer.View.CopyFromCPU(stream, ref MemoryMarshal.GetReference(powSpan), length); + orderBuffer.View.CopyFromCPU(stream, ref MemoryMarshal.GetReference(orderSpan), length); + resultBuffer.View.CopyFromCPU(stream, ref MemoryMarshal.GetReference(resultSpan), length); + statusBuffer.View.CopyFromCPU(stream, ref MemoryMarshal.GetReference(statusSpan), length); int pending; do { kernel(length, _divisorCyclesBatchSize, divisorBuffer.View, powBuffer.View, orderBuffer.View, resultBuffer.View, statusBuffer.View); - accelerator.Synchronize(); + stream.Synchronize(); - statusBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(statusSpan), length); + statusBuffer.View.CopyToCPU(stream, ref MemoryMarshal.GetReference(statusSpan), length); pending = 0; for (int i = 0; i < length; i++) { @@ -298,7 +299,7 @@ public sealed class DivisorCycleCache } while (pending > 0); - resultBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(resultSpan), length); + resultBuffer.View.CopyToCPU(stream, ref MemoryMarshal.GetReference(resultSpan), length); resultSpan.CopyTo(destination); if (rentedPow is not null) { diff --git a/PerfectNumbers.Core/Gpu/GpuContextPool.cs b/PerfectNumbers.Core/Gpu/GpuContextPool.cs index ae11bda..80e2adb 100644 --- a/PerfectNumbers.Core/Gpu/GpuContextPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuContextPool.cs @@ -4,6 +4,21 @@ using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu; +public sealed class GpuContextLease +{ + private GpuContextPool.PooledContext? _ctx; + + internal GpuContextLease(GpuContextPool.PooledContext ctx) + { + _ctx = ctx; + } + + public Context Context => _ctx!.Context; + public Accelerator Accelerator => _ctx!.Accelerator; + + public void Dispose() => GpuContextPool.Return(_ctx!); +} + public static class GpuContextPool { // The pool intentionally skips pre-loading ProcessEightBitWindows kernels so accelerator initialization @@ -14,12 +29,10 @@ public static class GpuContextPool public readonly Context Context; public readonly Accelerator Accelerator; - public PooledContext(KernelContainer kernels) + public PooledContext() { Context = SharedGpuContext.Context; Accelerator = SharedGpuContext.Accelerator; - GpuStaticTableInitializer.EnsureStaticTables(kernels, Accelerator); - kernels.Dispose(); // NOTE: Avoid loading/compiling any kernel here to prevent implicit // CL stream/queue creation during accelerator construction. // Some OpenCL drivers are fragile when a queue is created immediately @@ -29,11 +42,9 @@ public static class GpuContextPool public void Dispose() { // Ensure all cached GPU buffers for this accelerator are released. - NttGpuMath.ClearCaches(Accelerator); // Release PrimeTester GPU state for this accelerator (kernel/device primes) // Clear any per-accelerator cached resources to avoid releasing // after the accelerator is destroyed. - PrimeTester.ClearGpuCaches(Accelerator); // These resources are shared between GPU leases // Accelerator.Dispose(); @@ -53,7 +64,7 @@ public static class GpuContextPool } // Create a new accelerator when the pool does not have one available. - return new GpuContextLease(new PooledContext(GpuKernelPool.GetKernels())); + return new GpuContextLease(new PooledContext()); } @@ -71,8 +82,7 @@ public static class GpuContextPool PooledContext context; for (int i = 0; i < toCreate; i++) { - KernelContainer kernels = GpuKernelPool.GetKernels(); - context = new PooledContext(kernels); + context = new PooledContext(); contexts[i] = context; pool.Enqueue(context); } @@ -87,41 +97,12 @@ public static class GpuContextPool gpu.Dispose(); } - PrimeTester.DisposeGpuContexts(); WarmedGpuContextCount = 0; } - private static void Return(PooledContext ctx) + internal static void Return(PooledContext ctx) { - ctx.Accelerator.Synchronize(); + // ctx.Accelerator.Synchronize(); GpuPool.Enqueue(ctx); } - - public struct GpuContextLease - { - private PooledContext? _ctx; - private bool _disposed; - - internal GpuContextLease(PooledContext ctx) - { - _ctx = ctx; - _disposed = false; - } - - public Context Context => _ctx!.Context; - public Accelerator Accelerator => _ctx!.Accelerator; - - public void Dispose() - { - if (_disposed) - { - return; - } - - _disposed = true; - Return(_ctx!); - _ctx = null; - } - } -} - +} \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs index 1db24b1..fdbc1a5 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs @@ -1,9 +1,6 @@ -using System; using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; -using PerfectNumbers.Core; -using static PerfectNumbers.Core.Gpu.GpuContextPool; namespace PerfectNumbers.Core.Gpu; @@ -11,11 +8,13 @@ public sealed class GpuKernelLease { private static readonly ConcurrentQueue Pool = new(); - private GpuContextLease? _gpu; + private GpuContextLease _gpu; public KernelContainer Kernels; - private AcceleratorStream? _stream; + public AcceleratorStream Stream; - private GpuKernelLease(KernelContainer kernels) +#pragma warning disable CS8618 // _gpu and Stream will be set by Rent method + private GpuKernelLease(KernelContainer kernels) +#pragma warning restore CS8618 { Kernels = kernels; } @@ -28,35 +27,19 @@ public sealed class GpuKernelLease } lease._gpu = gpu; - lease.Kernels = kernels; - lease._stream = null; + lease.Kernels = kernels; + lease.Stream = gpu.Accelerator.CreateStream(); return lease; } public 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 Action, ArrayView> OrderKernel { get { var accel = Accelerator; // avoid capturing 'this' in lambda - return KernelContainer.InitOnce(ref Kernels!.Order, () => + return GpuKernelPool.InitOnce(ref Kernels!.Order, () => { var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernels.OrderKernelScan); var kernel = KernelUtil.GetKernel(loaded); @@ -70,7 +53,7 @@ public sealed class GpuKernelLease get { var accel = Accelerator; - return KernelContainer.InitOnce(ref Kernels!.Pow2Mod, () => + return GpuKernelPool.InitOnce(ref Kernels!.Pow2Mod, () => { var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernels.Pow2ModKernelScan); var kernel = KernelUtil.GetKernel(loaded); @@ -84,7 +67,7 @@ public sealed class GpuKernelLease get { var accel = Accelerator; - return KernelContainer.InitOnce(ref Kernels!.Incremental, () => + return GpuKernelPool.InitOnce(ref Kernels!.Incremental, () => { var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernels.IncrementalKernelScan); var kernel = KernelUtil.GetKernel(loaded); @@ -98,7 +81,7 @@ public sealed class GpuKernelLease get { var accel = Accelerator; - return KernelContainer.InitOnce(ref Kernels!.IncrementalOrder, () => + return GpuKernelPool.InitOnce(ref Kernels!.IncrementalOrder, () => { var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernels.IncrementalOrderKernelScan); var kernel = KernelUtil.GetKernel(loaded); @@ -112,7 +95,7 @@ public sealed class GpuKernelLease get { var accel = Accelerator; - return KernelContainer.InitOnce(ref Kernels!.Pow2ModOrder, () => + return GpuKernelPool.InitOnce(ref Kernels!.Pow2ModOrder, () => { var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModKernels.Pow2ModOrderKernelScan); var kernel = KernelUtil.GetKernel(loaded); @@ -121,16 +104,16 @@ public sealed class GpuKernelLease } } - public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D> SpecialMaxKernel + public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D> SpecialMaxKernel { get { var accel = Accelerator; - return KernelContainer.InitOnce(ref Kernels!.SpecialMax, () => + return GpuKernelPool.InitOnce(ref Kernels!.SpecialMax, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel); + var loaded = accel.LoadAutoGroupedStreamKernel, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel); var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>>(); + return kernel.CreateLauncherDelegate, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>>(); }); } } @@ -140,7 +123,7 @@ public sealed class GpuKernelLease get { var accel = Accelerator; - return KernelContainer.InitOnce(ref Kernels!.SmallPrimeFactor, () => + return GpuKernelPool.InitOnce(ref Kernels!.SmallPrimeFactor, () => { var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan); var kernel = KernelUtil.GetKernel(loaded); @@ -151,15 +134,10 @@ public sealed class GpuKernelLease public void Dispose() { - _stream?.Dispose(); - _stream = null; - - _gpu?.Dispose(); - _gpu = null; - - Kernels?.Dispose(); - GpuPrimeWorkLimiter.Release(); + Stream.Dispose(); + _gpu.Dispose(); Pool.Enqueue(this); + GpuPrimeWorkLimiter.Release(); } } diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 365b4d6..97fa8d8 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -1,4 +1,3 @@ -using PerfectNumbers.Core; using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; @@ -23,7 +22,8 @@ public readonly struct ResidueAutomatonArgs } } -public sealed class KernelContainer + +public sealed class KernelContainer : IDisposable { // Serializes first-time initialization of kernels/buffers per accelerator. public Action, ArrayView>? Order; @@ -38,42 +38,37 @@ public sealed class KernelContainer public Action, ArrayView1D>? Pow2ModOrder; public Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>? SmallPrimeFactor; - public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>? SpecialMax; + public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>? SpecialMax; // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length. public MemoryBuffer1D? SmallCycles; - public MemoryBuffer1D? SmallPrimeFactorsPrimes; - public MemoryBuffer1D? SmallPrimeFactorsSquares; - public MemoryBuffer1D? SmallPrimeFactorPrimeSlots; - public MemoryBuffer1D? SmallPrimeFactorExponentSlots; - public MemoryBuffer1D? SmallPrimeFactorCountSlot; - public MemoryBuffer1D? SmallPrimeFactorRemainingSlot; - public MemoryBuffer1D? SpecialMaxFactors; - public MemoryBuffer1D? SpecialMaxCandidates; - public MemoryBuffer1D? SpecialMaxResult; + public MemoryBuffer1D? SmallPrimesLastOne; - public MemoryBuffer1D? SmallPrimesPow2LastOne; public MemoryBuffer1D? SmallPrimesLastSeven; + public MemoryBuffer1D? SmallPrimesPow2LastOne; public MemoryBuffer1D? SmallPrimesPow2LastSeven; - public static T InitOnce(ref T? slot, Func factory) where T : class + public MemoryBuffer1D? SmallPrimeFactorsPrimes; + public MemoryBuffer1D? SmallPrimeFactorsSquares; + + public void Dispose() { - var current = Volatile.Read(ref slot); - if (current is not null) - { - return current; - } + // Order = null; + // Incremental = null; + // Pow2Mod = null; + // IncrementalOrder = null; + // Pow2ModOrder = null; + // SmallPrimeFactor = null; + // SpecialMax = null; - current = Volatile.Read(ref slot); - if (current is null) - { - current = factory(); - Volatile.Write(ref slot, current); - } - return current; + // SmallCycles?.Dispose(); + // SmallPrimesLastOne?.Dispose(); + // SmallPrimesLastSeven?.Dispose(); + // SmallPrimesPow2LastOne?.Dispose(); + // SmallPrimesPow2LastSeven?.Dispose(); + // SmallPrimeFactorsPrimes?.Dispose(); + // SmallPrimeFactorsSquares?.Dispose(); } - - public void Dispose() => GpuKernelPool.Return(this); } public readonly struct ResiduePrimeViews( @@ -101,49 +96,30 @@ public readonly struct SmallPrimeFactorTables( public ArrayView1D SquaresView => Squares.View; } -public readonly struct SmallPrimeFactorScratch( - MemoryBuffer1D primeSlots, - MemoryBuffer1D exponentSlots, - MemoryBuffer1D countSlot, - MemoryBuffer1D remainingSlot) -{ - public readonly MemoryBuffer1D PrimeSlots = primeSlots; - public readonly MemoryBuffer1D ExponentSlots = exponentSlots; - public readonly MemoryBuffer1D CountSlot = countSlot; - public readonly MemoryBuffer1D RemainingSlot = remainingSlot; - - public void Clear() - { - PrimeSlots.MemSetToZero(); - ExponentSlots.MemSetToZero(); - CountSlot.MemSetToZero(); - RemainingSlot.MemSetToZero(); - } -} - -public readonly struct SpecialMaxScratch( - MemoryBuffer1D factorValues, - MemoryBuffer1D candidateValues, - MemoryBuffer1D resultSlot) -{ - public readonly MemoryBuffer1D FactorValues = factorValues; - public readonly MemoryBuffer1D CandidateValues = candidateValues; - public readonly MemoryBuffer1D ResultSlot = resultSlot; - - public ArrayView1D FactorsView => FactorValues.View; - public ArrayView1D CandidatesView => CandidateValues.View; - public ArrayView1D ResultView => ResultSlot.View; -} - public class GpuKernelPool { - private static readonly ConcurrentQueue 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 readonly ConcurrentDictionary _kernels = new(); - public static KernelContainer GetKernels() => KernelCache.TryDequeue(out var container) - ? container - : new KernelContainer(); + public static T InitOnce(ref T? target, Func valueFactory) where T : class + { + if (target is { } value) + { + return value; + } - public static void Return(KernelContainer container) => KernelCache.Enqueue(container); + var newValue = valueFactory(); + return Interlocked.CompareExchange(ref target, newValue, null) ?? newValue; + } + + internal static KernelContainer GetKernels(Accelerator accelerator) + { + return _kernels.GetOrAdd(accelerator, static (accelerator) => + { + var kernels = new KernelContainer(); + PreloadStaticTables(kernels, accelerator); + return kernels; + }); + } internal static void PreloadStaticTables(KernelContainer kernels, Accelerator accelerator) { @@ -215,10 +191,10 @@ public class GpuKernelPool return new SmallPrimeFactorTables(primeBuffer, squareBuffer, (int)primeBuffer.Length); } - if (kernels.SmallPrimeFactorsPrimes is { } existingPrimes && kernels.SmallPrimeFactorsSquares is { } existingSquares) - { - return new SmallPrimeFactorTables(existingPrimes, existingSquares, (int)existingPrimes.Length); - } + // if (kernels.SmallPrimeFactorsPrimes is { } existingPrimes && kernels.SmallPrimeFactorsSquares is { } existingSquares) + // { + // return new SmallPrimeFactorTables(existingPrimes, existingSquares, (int)existingPrimes.Length); + // } var hostPrimes = PrimesGenerator.SmallPrimes; var hostSquares = PrimesGenerator.SmallPrimesPow2; @@ -234,74 +210,15 @@ public class GpuKernelPool return new SmallPrimeFactorTables(devicePrimes, deviceSquares, hostPrimes.Length); } - public static SmallPrimeFactorScratch EnsureSmallPrimeFactorScratch(KernelContainer kernels, Accelerator accelerator, int slotCount) - { - MemoryBuffer1D? primeSlots = kernels.SmallPrimeFactorPrimeSlots; - if (primeSlots is null || primeSlots.Length < slotCount) - { - primeSlots?.Dispose(); - primeSlots = accelerator.Allocate1D(slotCount); - kernels.SmallPrimeFactorPrimeSlots = primeSlots; - } - MemoryBuffer1D? exponentSlots = kernels.SmallPrimeFactorExponentSlots; - if (exponentSlots is null || exponentSlots.Length < slotCount) - { - exponentSlots?.Dispose(); - exponentSlots = accelerator.Allocate1D(slotCount); - kernels.SmallPrimeFactorExponentSlots = exponentSlots; - } - MemoryBuffer1D? countSlot = kernels.SmallPrimeFactorCountSlot; - if (countSlot is null) - { - countSlot = accelerator.Allocate1D(1); - kernels.SmallPrimeFactorCountSlot = countSlot; - } - MemoryBuffer1D? remainingSlot = kernels.SmallPrimeFactorRemainingSlot; - if (remainingSlot is null) - { - remainingSlot = accelerator.Allocate1D(1); - kernels.SmallPrimeFactorRemainingSlot = remainingSlot; - } - return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!); - } - - public static SpecialMaxScratch EnsureSpecialMaxScratch(KernelContainer kernels, Accelerator accelerator, int factorCapacity) - { - MemoryBuffer1D? factorValues = kernels.SpecialMaxFactors; - if (factorValues is null || factorValues.Length < factorCapacity) - { - factorValues?.Dispose(); - factorValues = accelerator.Allocate1D(factorCapacity); - kernels.SpecialMaxFactors = factorValues; - } - - MemoryBuffer1D? candidateValues = kernels.SpecialMaxCandidates; - if (candidateValues is null || candidateValues.Length < factorCapacity) - { - candidateValues?.Dispose(); - candidateValues = accelerator.Allocate1D(factorCapacity); - kernels.SpecialMaxCandidates = candidateValues; - } - - MemoryBuffer1D? resultSlot = kernels.SpecialMaxResult; - if (resultSlot is null) - { - resultSlot = accelerator.Allocate1D(1); - kernels.SpecialMaxResult = resultSlot; - } - - return new SpecialMaxScratch(factorValues!, candidateValues!, resultSlot!); - } - - public static GpuKernelLease GetKernel(bool useGpuOrder) + public static GpuKernelLease GetKernel() { GpuPrimeWorkLimiter.Acquire(); var gpu = Rent(); - var kernels = GetKernels(); + var kernels = GetKernels(gpu.Accelerator); return GpuKernelLease.Rent(gpu, kernels); } @@ -311,12 +228,10 @@ public class GpuKernelPool /// Action to run with (Accelerator, Stream). public static void Run(Action action) { - var lease = GetKernel(useGpuOrder: true); + var lease = GetKernel(); var accelerator = lease.Accelerator; - var stream = accelerator.CreateStream(); + var stream = lease.Stream; action(accelerator, stream); - stream.Dispose(); lease.Dispose(); - } } diff --git a/PerfectNumbers.Core/Gpu/GpuScratchBufferPool.cs b/PerfectNumbers.Core/Gpu/GpuScratchBufferPool.cs new file mode 100644 index 0000000..759b2a7 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/GpuScratchBufferPool.cs @@ -0,0 +1,134 @@ +using System.Collections.Concurrent; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +public sealed class ScratchBuffer +{ + public MemoryBuffer1D SmallPrimeFactorPrimeSlots { get; private set; } + public MemoryBuffer1D SmallPrimeFactorExponentSlots { get; private set; } + public MemoryBuffer1D SmallPrimeFactorCountSlot { get; private set; } + public MemoryBuffer1D SmallPrimeFactorRemainingSlot { get; private set; } + + public MemoryBuffer1D SpecialMaxFactors { get; private set; } + public MemoryBuffer1D SpecialMaxCandidates { get; private set; } + public MemoryBuffer1D SpecialMaxResult { get; private set; } + + internal readonly Accelerator _accelerator; + + public ScratchBuffer(Accelerator accelerator, int smallPrimeFactorSlotCount, int specialMaxFactorCapacity) + { + _accelerator = accelerator; + + // lock (accelerator) + // { + SmallPrimeFactorPrimeSlots = accelerator.Allocate1D(smallPrimeFactorSlotCount); + SmallPrimeFactorExponentSlots = accelerator.Allocate1D(smallPrimeFactorSlotCount); + SmallPrimeFactorCountSlot = accelerator.Allocate1D(1); + SmallPrimeFactorRemainingSlot = accelerator.Allocate1D(1); + + SpecialMaxFactors = accelerator.Allocate1D(specialMaxFactorCapacity); + SpecialMaxCandidates = accelerator.Allocate1D(specialMaxFactorCapacity); + SpecialMaxResult = accelerator.Allocate1D(1); + // } + } + + public void ResizeSmallPrimeFactorSlots(int newSize) + { + if (SmallPrimeFactorPrimeSlots.Length < newSize) + { + SmallPrimeFactorPrimeSlots.Dispose(); + SmallPrimeFactorExponentSlots.Dispose(); + + SmallPrimeFactorPrimeSlots = _accelerator.Allocate1D(newSize); + SmallPrimeFactorExponentSlots = _accelerator.Allocate1D(newSize); + } + // else + // { + // SmallPrimeFactorPrimeSlots.MemSetToZero(); + // SmallPrimeFactorExponentSlots.MemSetToZero(); + // } + } + + public void ResizeSpecialMaxFactors(int newSize) + { + if (SpecialMaxFactors.Length < newSize) + { + SpecialMaxFactors.Dispose(); + SpecialMaxCandidates.Dispose(); + + SpecialMaxFactors = _accelerator.Allocate1D(newSize); + SpecialMaxCandidates = _accelerator.Allocate1D(newSize); + } + // else + // { + // SpecialMaxFactors.MemSetToZero(); + // SpecialMaxCandidates.MemSetToZero(); + // } + } +} + +public static class GpuScratchBufferPool +{ + [ThreadStatic] + private static Dictionary>? _pools; + // private static readonly ConcurrentDictionary> _pools = new(); + private const int DefaultSmallPrimeFactorSlotCount = 64; // From PrimeOrderCalculator.Gpu.cs + private const int DefaultSpecialMaxFactorCapacity = 1024; // A reasonable default, will be resized if needed + + public static ScratchBuffer Rent(Accelerator accelerator, int smallPrimeFactorSlotCount, int specialMaxFactorCapacity) + { + const int WarmUpCapacity = 100; + // var pool = _pools.GetOrAdd(accelerator, _ => new ConcurrentQueue()); + if (_pools is not { }) + { + _pools = new(WarmUpCapacity); + WarmUp(WarmUpCapacity, DefaultSmallPrimeFactorSlotCount, DefaultSpecialMaxFactorCapacity); + } + + var pools = _pools; + if (!pools.TryGetValue(accelerator, out var pool)) + { + pool = new(); + pools[accelerator] = pool; + } + + // var pool = _pools.GetOrAdd(accelerator, _ => new ConcurrentQueue()); + if (pool.TryDequeue(out var buffer)) + { + if (buffer.SmallPrimeFactorPrimeSlots.Length < smallPrimeFactorSlotCount || buffer.SpecialMaxFactors.Length < specialMaxFactorCapacity) + { + throw new InvalidOperationException(); + // lock (accelerator) + // { + // Console.WriteLine($"Resizing GPU scratch buffer from pool ({buffer.SmallPrimeFactorPrimeSlots.Length} / {smallPrimeFactorSlotCount}), ({buffer.SpecialMaxFactors.Length}/{specialMaxFactorCapacity})"); + // buffer.ResizeSmallPrimeFactorSlots(smallPrimeFactorSlotCount); + // buffer.ResizeSpecialMaxFactors(specialMaxFactorCapacity); + // } + } + return buffer; + } + + return new ScratchBuffer(accelerator, smallPrimeFactorSlotCount, specialMaxFactorCapacity); + } + + public static void Return(ScratchBuffer buffer) + { + if (_pools!.TryGetValue(buffer._accelerator, out var pool)) + { + pool.Enqueue(buffer); + } + } + + public static void WarmUp(int count, int smallPrimeFactorSlots, int specialMaxFactorCapacity) + { + Accelerator accelerator = SharedGpuContext.Accelerator; + + for (int i = 0; i < count; i++) + { + var buffer = Rent(accelerator, smallPrimeFactorSlots, specialMaxFactorCapacity); + Return(buffer); + } + } +} \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs index 3b71ec2..6ea5dd0 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs @@ -1,6 +1,5 @@ using ILGPU; using ILGPU.Runtime; -using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -13,7 +12,7 @@ internal static partial class PrimeOrderGpuHeuristics int factorCount, MontgomeryDivisorData divisor, ArrayView1D candidates, - ArrayView1D resultOut) + ArrayView1D resultOut) { if (index != 0) { @@ -62,7 +61,7 @@ internal static partial class PrimeOrderGpuHeuristics var stepper = new ExponentRemainderStepperGpu(divisorPartial); if (stepper.InitializeIsUnityGpu(candidates[0])) { - resultOut[0] = 0; + Atomic.Exchange(ref resultOut[0], 0); return; } @@ -70,11 +69,11 @@ internal static partial class PrimeOrderGpuHeuristics { if (stepper.ComputeNextIsUnityGpu(candidates[i])) { - resultOut[0] = 0; + Atomic.Exchange(ref resultOut[0], 0); return; } } - resultOut[0] = 1; + Atomic.Exchange(ref resultOut[0], 1); } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 19b08c9..b4988a0 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -18,19 +18,19 @@ internal static class PrimeTesterKernels { // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise. public static void SmallPrimeSieveKernel( - Index1D index, - ArrayView numbers, - ArrayView smallPrimes, - ArrayView smallPrimesLastOne, - ArrayView smallPrimesLastSeven, - ArrayView smallPrimesLastThree, - ArrayView smallPrimesLastNine, - ArrayView smallPrimesPow2, - ArrayView smallPrimesPow2LastOne, - ArrayView smallPrimesPow2LastSeven, - ArrayView smallPrimesPow2LastThree, - ArrayView smallPrimesPow2LastNine, - ArrayView results) + Index1D index, + ArrayView numbers, + ArrayView smallPrimes, + ArrayView smallPrimesLastOne, + ArrayView smallPrimesLastSeven, + ArrayView smallPrimesLastThree, + ArrayView smallPrimesLastNine, + ArrayView smallPrimesPow2, + ArrayView smallPrimesPow2LastOne, + ArrayView smallPrimesPow2LastSeven, + ArrayView smallPrimesPow2LastThree, + ArrayView smallPrimesPow2LastNine, + ArrayView results) { ulong n = numbers[index]; byte result = 1; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 25391b0..0648c2c 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,3 +1,4 @@ +using PerfectNumbers.Core.Gpu; using System; using System.Buffers; using System.Collections.Concurrent; @@ -25,7 +26,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>> _kernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _resourcePool = new(); - private readonly ConcurrentBag _acceleratorPool = new(); + private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); private Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => @@ -476,7 +477,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession { private readonly MersenneNumberDivisorByDivisorGpuTester _owner; - private readonly GpuContextPool.GpuContextLease _lease; + private readonly GpuContextLease _lease; private readonly Accelerator _accelerator; private Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> _kernel = null!; private MemoryBuffer1D _divisorBuffer = null!; @@ -650,7 +651,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } - private GpuContextPool.GpuContextLease RentAccelerator() + private GpuContextLease RentAccelerator() { if (_acceleratorPool.TryTake(out var lease)) { @@ -660,7 +661,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return GpuContextPool.Rent(); } - private void ReturnAccelerator(GpuContextPool.GpuContextLease lease) => _acceleratorPool.Add(lease); + private void ReturnAccelerator(GpuContextLease lease) => _acceleratorPool.Add(lease); private BatchResources RentBatchResources(Accelerator accelerator, int capacity) { @@ -899,8 +900,4 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() => CreateDivisorSession(); -} - - - - +} \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index eaffd13..3114bcb 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -38,19 +38,21 @@ public sealed class MersenneNumberDivisorGpuTester public bool IsDivisible(ulong exponent, in ReadOnlyGpuUInt128 divisor) { var gpu = GpuContextPool.Rent(); - var accelerator = gpu.Accelerator; + var accelerator = gpu.Accelerator; + var stream = accelerator.CreateStream(); var kernel = GetKernel(accelerator); var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D(1)); // There is no point in clearing this buffer. We always override item [0] and never use it beyond item [0] // resultBuffer.MemSetToZero(); kernel(1, exponent, divisor, resultBuffer.View); - accelerator.Synchronize(); + stream.Synchronize(); Span result = stackalloc byte[1]; - resultBuffer.View.CopyToCPU(ref result[0], 1); + resultBuffer.View.CopyToCPU(stream, ref result[0], 1); bool divisible = result[0] != 0; result[0] = 0; - resultBuffer.View.CopyFromCPU(ref result[0], 1); - accelerator.Synchronize(); + resultBuffer.View.CopyFromCPU(stream, ref result[0], 1); + stream.Synchronize(); + stream.Dispose(); gpu.Dispose(); return divisible; } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs index a9467dc..76c8b72 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs @@ -16,7 +16,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u { throw new NotImplementedException($"GPU incremental scanning requires the device cycle heuristics implementation (kernel: {_kernelType}, GPU order: {_useGpuOrder})."); - var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); + var gpuLease = GpuKernelPool.GetKernel(); var accelerator = gpuLease.Accelerator; var stream = gpuLease.Stream; int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs index 8c72652..89d86f8 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs @@ -9,7 +9,7 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { - var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); + var gpuLease = GpuKernelPool.GetKernel(); var accelerator = gpuLease.Accelerator; var stream = gpuLease.Stream; int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy and avoids TDR diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index fdd1e84..5c36843 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -1,3 +1,4 @@ +using PerfectNumbers.Core.Gpu; using System; using System.Buffers; using System.Collections.Concurrent; @@ -13,7 +14,7 @@ namespace PerfectNumbers.Core.Gpu; public class MersenneNumberResidueGpuTester(bool useGpuOrder) { private readonly bool _useGpuOrder = useGpuOrder; - private readonly ConcurrentBag _acceleratorPool = new(); + private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentDictionary, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>> _pow2ModKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); @@ -35,7 +36,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) var stream = accelerator.CreateStream(); var kernel = GetPow2ModKernel(accelerator); var resources = RentResources(accelerator, GpuConstants.ScanBatchSize); - var gpuKernels = GpuKernelPool.GetKernels(); + var gpuKernels = GpuKernelPool.GetKernels(accelerator); var orderBuffer = resources.OrderBuffer; ulong[] orderArray = resources.OrderArray; @@ -69,7 +70,8 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL, kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2); - stream.Synchronize(); + stream.Synchronize(); + orderBuffer.View.CopyToCPU(ref ordersRef, currentSize); if (!Volatile.Read(ref isPrime)) { @@ -127,7 +129,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) loaded.Method); }); - private GpuContextPool.GpuContextLease RentAccelerator() + private GpuContextLease RentAccelerator() { if (_acceleratorPool.TryTake(out var lease)) { @@ -137,11 +139,11 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) return GpuContextPool.Rent(); } - private void ReturnAccelerator(GpuContextPool.GpuContextLease lease) => _acceleratorPool.Add(lease); + private void ReturnAccelerator(GpuContextLease lease) => _acceleratorPool.Add(lease); private ResidueResources RentResources(Accelerator accelerator, int capacity) { - var bag = _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()); + var bag = _resourcePools.GetOrAdd(accelerator, static _ => []); while (bag.TryTake(out var resources)) { if (resources.Capacity >= capacity) @@ -195,4 +197,4 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) } -} +} \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 8664d8b..931bbb9 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -198,7 +198,7 @@ internal static partial class PrimeOrderGpuHeuristics try { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var lease = GpuKernelPool.GetKernel(); Accelerator accelerator = lease.Accelerator; AcceleratorStream stream = lease.Stream; @@ -364,7 +364,7 @@ internal static partial class PrimeOrderGpuHeuristics { order = 0UL; - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var lease = GpuKernelPool.GetKernel(); Accelerator accelerator = lease.Accelerator; AcceleratorStream stream = lease.Stream; @@ -526,7 +526,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool TryComputeOnGpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var lease = GpuKernelPool.GetKernel(); Accelerator accelerator = lease.Accelerator; AcceleratorStream stream = lease.Stream; var kernel = GetPow2ModKernel(accelerator); @@ -581,7 +581,7 @@ internal static partial class PrimeOrderGpuHeuristics GpuUInt128[]? rentedExponents = null; GpuUInt128[]? rentedResults = null; - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var lease = GpuKernelPool.GetKernel(); try { diff --git a/PerfectNumbers.Core/Gpu/SharedGpuContext.cs b/PerfectNumbers.Core/Gpu/SharedGpuContext.cs index ec5c4cb..b594d66 100644 --- a/PerfectNumbers.Core/Gpu/SharedGpuContext.cs +++ b/PerfectNumbers.Core/Gpu/SharedGpuContext.cs @@ -5,7 +5,7 @@ using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu; -internal static class SharedGpuContext +public static class SharedGpuContext { private static readonly Lazy s_context = new(() => ILGPU.Context.CreateDefault(), LazyThreadSafetyMode.ExecutionAndPublication); private static readonly Lazy s_device = new(() => s_context.Value.GetPreferredDevice(false), LazyThreadSafetyMode.ExecutionAndPublication); diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 4389732..7194f02 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -2,6 +2,7 @@ using Open.Numeric.Primes; using ILGPU.Runtime; using System.Runtime.CompilerServices; using PerfectNumbers.Core.Gpu; +using System.Runtime.InteropServices; namespace PerfectNumbers.Core; @@ -227,6 +228,7 @@ public sealed class HeuristicCombinedPrimeTester GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(1); var accelerator = gpu.Accelerator; + var stream = gpu.Stream; var kernel = PrimeTester.PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.HeuristicTrialDivisionKernel; var flagView1D = gpu.HeuristicFlag.View; var flagView = flagView1D.AsContiguous(); @@ -243,8 +245,9 @@ public sealed class HeuristicCombinedPrimeTester maxDivisorSquare, HeuristicGpuDivisorTableKind.Combined, gpu.HeuristicGpuTables); - accelerator.Synchronize(); - flagView1D.CopyToCPU(ref compositeFlag, 1); + stream.Synchronize(); + + flagView1D.CopyToCPU(stream, ref compositeFlag, 1); compositeDetected = compositeFlag != 0; gpu.Dispose(); diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index 67c97b1..18ca587 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -243,6 +243,7 @@ public sealed class HeuristicPrimeTester GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(1); var accelerator = gpu.Accelerator; + var stream = gpu.Stream; var kernel = PrimeTester.PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.HeuristicTrialDivisionKernel; var flagView1D = gpu.HeuristicFlag.View; var flagView = flagView1D.AsContiguous(); @@ -260,8 +261,9 @@ public sealed class HeuristicPrimeTester maxDivisorSquare, HeuristicGpuDivisorTableKind.GroupA, gpu.HeuristicGpuTables); - accelerator.Synchronize(); - flagView1D.CopyToCPU(ref compositeFlag, 1); + stream.Synchronize(); + + flagView1D.CopyToCPU(stream, ref compositeFlag, 1); compositeDetected = compositeFlag != 0; if (!compositeDetected) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 9e7b026..0662b5e 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -13,364 +13,364 @@ namespace PerfectNumbers.Core; public class MersenneDivisorCycles { - private List<(ulong divisor, ulong cycleLength)> _table = []; - // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown - private ulong[]? _smallCycles; - public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles(); + private List<(ulong divisor, ulong cycleLength)> _table = []; + // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown + private ulong[]? _smallCycles; + public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles(); - public MersenneDivisorCycles() - { - } + public MersenneDivisorCycles() + { + } - public void LoadFrom(string path) - { - EnsureSmallBuffer(); - List<(ulong divisor, ulong cycle)> cycles = []; - using Stream outputStream = new FileStream(path, FileMode.Open, FileAccess.Read, FileShare.Read, BufferSize10M, useAsync: false); - foreach (var (d, c) in EnumerateStream(outputStream)) - { - cycles.Add((d, c)); - if (d <= PerfectNumberConstants.MaxQForDivisorCycles) - { - _smallCycles![(int)d] = c == 0UL ? 1UL : c; - } - } + public void LoadFrom(string path) + { + EnsureSmallBuffer(); + List<(ulong divisor, ulong cycle)> cycles = []; + using Stream outputStream = new FileStream(path, FileMode.Open, FileAccess.Read, FileShare.Read, BufferSize10M, useAsync: false); + foreach (var (d, c) in EnumerateStream(outputStream)) + { + cycles.Add((d, c)); + if (d <= PerfectNumberConstants.MaxQForDivisorCycles) + { + _smallCycles![(int)d] = c == 0UL ? 1UL : c; + } + } - cycles.Sort((a, b) => a.divisor < b.divisor ? -1 : (a.divisor > b.divisor ? 1 : 0)); - _table = cycles; - } + cycles.Sort((a, b) => a.divisor < b.divisor ? -1 : (a.divisor > b.divisor ? 1 : 0)); + _table = cycles; + } - public const int SmallDivisorsMax = PerfectNumberConstants.MaxQForDivisorCycles; + public const int SmallDivisorsMax = PerfectNumberConstants.MaxQForDivisorCycles; - // Provides a snapshot array of small cycles [0..PerfectNumberConstants.MaxQForDivisorCycles]. Index is the divisor. - // The returned array is a copy safe for use across threads and device uploads. - public ulong[] ExportSmallCyclesSnapshot() - { - EnsureSmallBuffer(); - ulong[] snapshot = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; - Array.Copy(_smallCycles!, snapshot, snapshot.Length); - return snapshot; - } + // Provides a snapshot array of small cycles [0..PerfectNumberConstants.MaxQForDivisorCycles]. Index is the divisor. + // The returned array is a copy safe for use across threads and device uploads. + public ulong[] ExportSmallCyclesSnapshot() + { + EnsureSmallBuffer(); + ulong[] snapshot = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + Array.Copy(_smallCycles!, snapshot, snapshot.Length); + return snapshot; + } - public static bool CycleEqualsExponent(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) - { - if (CycleEqualsExponentForMersenneCandidate(divisor, divisorData, exponent)) - { - return true; - } + public static bool CycleEqualsExponent(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) + { + if (CycleEqualsExponentForMersenneCandidate(divisor, divisorData, exponent)) + { + return true; + } - if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) - { - return false; - } + if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + { + return false; + } - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - var small = Shared._smallCycles; - if (small is not null) - { - ulong cached = small[(int)divisor]; - if (cached != 0UL) - { - return cached == exponent; - } - } - } + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + var small = Shared._smallCycles; + if (small is not null) + { + ulong cached = small[(int)divisor]; + if (cached != 0UL) + { + return cached == exponent; + } + } + } - if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL) - { - return heuristicCycle == exponent; - } + if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL) + { + return heuristicCycle == exponent; + } - if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL) - { - return cycleLength == exponent; - } + if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL) + { + return cycleLength == exponent; + } - return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; - } + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; + } - public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) - { - // EvenPerfectBitScanner only calls this helper with divisors generated from 2 * k * p + 1 and primes greater than one, - // so the validation guard would never trigger in production. Keep it commented out to leave the hot path branch-free. - // if (!IsValidMersenneDivisorCandidate(divisor, exponent)) - // { - // return false; - // } + public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) + { + // EvenPerfectBitScanner only calls this helper with divisors generated from 2 * k * p + 1 and primes greater than one, + // so the validation guard would never trigger in production. Keep it commented out to leave the hot path branch-free. + // if (!IsValidMersenneDivisorCandidate(divisor, exponent)) + // { + // return false; + // } - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - MersenneDivisorCycles shared = Shared; - ulong[]? small = shared._smallCycles; - if (small is not null) - { - ulong cached = small[(int)divisor]; - if (cached != 0UL) - { - return cached == exponent; - } - } - } + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + MersenneDivisorCycles shared = Shared; + ulong[]? small = shared._smallCycles; + if (small is not null) + { + ulong cached = small[(int)divisor]; + if (cached != 0UL) + { + return cached == exponent; + } + } + } - return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; - } + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; + } - private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) - { - if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) - { - return false; - } + private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) + { + if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + { + return false; + } - UInt128 step = (UInt128)exponent << 1; - if (step == UInt128.Zero) - { - return false; - } + UInt128 step = (UInt128)exponent << 1; + if (step == UInt128.Zero) + { + return false; + } - UInt128 adjusted = (UInt128)divisor - UInt128.One; - return adjusted % step == UInt128.Zero; - } + UInt128 adjusted = (UInt128)divisor - UInt128.One; + return adjusted % step == UInt128.Zero; + } - public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor) - { - // Binary pairs: (ulong divisor, ulong cycle) - using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true); - Stream stream = reader.BaseStream; - byte[] buffer = new byte[16]; - while (true) - { - int read = 0; - Span span = buffer; - while (read < span.Length) - { - int chunk = stream.Read(span[read..]); - if (chunk == 0) - { - if (read == 0) - { - yield break; - } + public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor) + { + // Binary pairs: (ulong divisor, ulong cycle) + using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true); + Stream stream = reader.BaseStream; + byte[] buffer = new byte[16]; + while (true) + { + int read = 0; + Span span = buffer; + while (read < span.Length) + { + int chunk = stream.Read(span[read..]); + if (chunk == 0) + { + if (read == 0) + { + yield break; + } - throw new EndOfStreamException("Incomplete divisor-cycle entry."); - } + throw new EndOfStreamException("Incomplete divisor-cycle entry."); + } - read += chunk; - } + read += chunk; + } - ulong divisor = BinaryPrimitives.ReadUInt64LittleEndian(span[..8]); - ulong cycle = BinaryPrimitives.ReadUInt64LittleEndian(span[8..]); - yield return (divisor, cycle); - } - } + ulong divisor = BinaryPrimitives.ReadUInt64LittleEndian(span[..8]); + ulong cycle = BinaryPrimitives.ReadUInt64LittleEndian(span[8..]); + yield return (divisor, cycle); + } + } - public ulong GetCycle(ulong divisor, in MontgomeryDivisorData divisorData) - { - // Fast-path: in-memory array for small divisors - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - var arr = _smallCycles; - if (arr is not null) - { - ulong cached = arr[(int)divisor]; - if (cached != 0UL) - { - return cached; - } - } - } + public ulong GetCycle(ulong divisor, in MontgomeryDivisorData divisorData) + { + // Fast-path: in-memory array for small divisors + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + var arr = _smallCycles; + if (arr is not null) + { + ulong cached = arr[(int)divisor]; + if (cached != 0UL) + { + return cached; + } + } + } - // Binary search for divisor in sorted _table - int mid, right, left = 0; - ulong d; - List<(ulong divisor, ulong cycleLength)> cycles = _table; + // Binary search for divisor in sorted _table + int mid, right, left = 0; + ulong d; + List<(ulong divisor, ulong cycleLength)> cycles = _table; - right = cycles.Count - 1; - while (left <= right) - { - mid = left + ((right - left) >> 1); - d = cycles[mid].divisor; + right = cycles.Count - 1; + while (left <= right) + { + mid = left + ((right - left) >> 1); + d = cycles[mid].divisor; - if (d == divisor) - { - return cycles[mid].cycleLength; - } + if (d == divisor) + { + return cycles[mid].cycleLength; + } - if (d < divisor) - { - left = mid + 1; - } - else - { - right = mid - 1; - } - } + if (d < divisor) + { + left = mid + 1; + } + else + { + right = mid - 1; + } + } - // TODO: Route this miss to an ephemeral single-cycle computation on the device selected by the current - // configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot - // so we respect the large-p memory limits and avoid extra cache locks. - return CalculateCycleLength(divisor, divisorData); - } + // TODO: Route this miss to an ephemeral single-cycle computation on the device selected by the current + // configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot + // so we respect the large-p memory limits and avoid extra cache locks. + return CalculateCycleLength(divisor, divisorData); + } - public static UInt128 GetCycle(UInt128 divisor) - { - // For divisor = 2^k, cycle is 1 - if ((divisor & (divisor - UInt128.One)) == UInt128.Zero) - { - return UInt128.One; - } + public static UInt128 GetCycle(UInt128 divisor) + { + // For divisor = 2^k, cycle is 1 + if ((divisor & (divisor - UInt128.One)) == UInt128.Zero) + { + return UInt128.One; + } - if (divisor <= ulong.MaxValue) - { - ulong divisor64 = (ulong)divisor; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64); - ulong cycle = CalculateCycleLength(divisor64, divisorData); - return (UInt128)cycle; - } + if (divisor <= ulong.MaxValue) + { + ulong divisor64 = (ulong)divisor; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64); + ulong cycle = CalculateCycleLength(divisor64, divisorData); + return (UInt128)cycle; + } - UInt128 wideOrder = PrimeOrderCalculator.Calculate( - divisor, - previousOrder: null, - PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, - PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); - if (wideOrder != UInt128.Zero) - { - return wideOrder; - } + UInt128 wideOrder = PrimeOrderCalculator.Calculate( + divisor, + previousOrder: null, + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); + if (wideOrder != UInt128.Zero) + { + return wideOrder; + } - // Otherwise, find order of 2 mod divisor - // TODO: Port this UInt128 path to the unrolled-hex cycle calculator so wide - // divisors stop relying on the slow shift-and-subtract loop measured in the - // GPU benchmarks. - UInt128 order = UInt128.One; - UInt128 pow = UInt128Numbers.Two; - while (pow != UInt128.One) - { - pow <<= 1; - if (pow >= divisor) - { - pow -= divisor; - } + // Otherwise, find order of 2 mod divisor + // TODO: Port this UInt128 path to the unrolled-hex cycle calculator so wide + // divisors stop relying on the slow shift-and-subtract loop measured in the + // GPU benchmarks. + UInt128 order = UInt128.One; + UInt128 pow = UInt128Numbers.Two; + while (pow != UInt128.One) + { + pow <<= 1; + if (pow >= divisor) + { + pow -= divisor; + } - order++; - } + order++; + } - return order; - } + return order; + } - public static bool TryCalculateCycleLengthForExponentCpu( - ulong divisor, - ulong exponent, - in MontgomeryDivisorData divisorData, - out ulong cycleLength, - out bool primeOrderFailed) - { - cycleLength = 0UL; - primeOrderFailed = false; + public static bool TryCalculateCycleLengthForExponentCpu( + ulong divisor, + ulong exponent, + in MontgomeryDivisorData divisorData, + out ulong cycleLength, + out bool primeOrderFailed) + { + cycleLength = 0UL; + primeOrderFailed = false; - // EvenPerfectBitScanner never feeds powers of two into this helper, so leave the branch commented out to avoid redundant work. - // if ((divisor & (divisor - 1UL)) == 0UL) - // { - // cycleLength = 1UL; - // return true; - // } + // EvenPerfectBitScanner never feeds powers of two into this helper, so leave the branch commented out to avoid redundant work. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // cycleLength = 1UL; + // return true; + // } - // The scanner only calls this path with prime exponents greater than one and odd divisors above three, keeping the small-parameter guard unreachable. - // if (divisor <= 3UL || exponent <= 1UL) - // { - // cycleLength = CalculateCycleLength(divisor, divisorData); - // return true; - // } + // The scanner only calls this path with prime exponents greater than one and odd divisors above three, keeping the small-parameter guard unreachable. + // if (divisor <= 3UL || exponent <= 1UL) + // { + // cycleLength = CalculateCycleLength(divisor, divisorData); + // return true; + // } - ulong computedOrder = PrimeOrderCalculator.Calculate( - divisor, - previousOrder: null, - divisorData, - PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, - PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); - if (computedOrder != 0UL) - { - cycleLength = computedOrder; - return true; - } + ulong computedOrder = PrimeOrderCalculator.Calculate( + divisor, + previousOrder: null, + divisorData, + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); + if (computedOrder != 0UL) + { + cycleLength = computedOrder; + return true; + } - primeOrderFailed = true; + primeOrderFailed = true; - ulong phi = divisor - 1UL; - // A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without - // branching in the hot loop. - // if (phi == 0UL) - // { - // return false; - // } + ulong phi = divisor - 1UL; + // A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without + // branching in the hot loop. + // if (phi == 0UL) + // { + // return false; + // } - int twoCount = BitOperations.TrailingZeroCount(phi); - ulong reducedPhi = phi >> twoCount; - if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL) - { - return false; - } + int twoCount = BitOperations.TrailingZeroCount(phi); + ulong reducedPhi = phi >> twoCount; + if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL) + { + return false; + } - ulong k = reducedPhi / exponent; - // k == 0 would indicate reducedPhi < exponent, which the EvenPerfectBitScanner call flow prevents. Document the - // assumption instead of branching. - // if (k == 0UL) - // { - // return false; - // } + ulong k = reducedPhi / exponent; + // k == 0 would indicate reducedPhi < exponent, which the EvenPerfectBitScanner call flow prevents. Document the + // assumption instead of branching. + // if (k == 0UL) + // { + // return false; + // } - Dictionary factorCounts = ThreadStaticPools.RentFactorCountDictionary(); - bool success = false; + Dictionary factorCounts = ThreadStaticPools.RentFactorCountDictionary(); + bool success = false; - if (twoCount > 0) - { - factorCounts[2UL] = twoCount; - } + if (twoCount > 0) + { + factorCounts[2UL] = twoCount; + } - bool exponentAccumulated = AccumulateFactors(exponent, factorCounts); - bool factorsAccumulated = exponentAccumulated && AccumulateFactors(k, factorCounts); + bool exponentAccumulated = AccumulateFactors(exponent, factorCounts); + bool factorsAccumulated = exponentAccumulated && AccumulateFactors(k, factorCounts); - if (factorsAccumulated) - { - cycleLength = ReduceOrder(divisorData, phi, factorCounts); - success = true; - primeOrderFailed = false; - } + if (factorsAccumulated) + { + cycleLength = ReduceOrder(divisorData, phi, factorCounts); + success = true; + primeOrderFailed = false; + } - factorCounts.Clear(); - ThreadStaticPools.ReturnFactorCountDictionary(factorCounts); - return success; - } + factorCounts.Clear(); + ThreadStaticPools.ReturnFactorCountDictionary(factorCounts); + return success; + } - private static bool AccumulateFactors( - ulong value, - Dictionary counts) - { - if (value <= 1UL) - { - return true; - } + private static bool AccumulateFactors( + ulong value, + Dictionary counts) + { + if (value <= 1UL) + { + return true; + } - Dictionary scratch = ThreadStaticPools.RentFactorScratchDictionary(); - bool success = TryFactorIntoCountsInternal(value, scratch); - if (!success) - { - scratch.Clear(); - ThreadStaticPools.ReturnFactorScratchDictionary(scratch); - return false; - } + Dictionary scratch = ThreadStaticPools.RentFactorScratchDictionary(); + bool success = TryFactorIntoCountsInternal(value, scratch); + if (!success) + { + scratch.Clear(); + ThreadStaticPools.ReturnFactorScratchDictionary(scratch); + return false; + } - for (Dictionary.Enumerator enumerator = scratch.GetEnumerator(); enumerator.MoveNext();) - { - KeyValuePair pair = enumerator.Current; - AddFactor(counts, pair.Key, pair.Value); - } + for (Dictionary.Enumerator enumerator = scratch.GetEnumerator(); enumerator.MoveNext();) + { + KeyValuePair pair = enumerator.Current; + AddFactor(counts, pair.Key, pair.Value); + } - scratch.Clear(); - ThreadStaticPools.ReturnFactorScratchDictionary(scratch); - return true; - } + scratch.Clear(); + ThreadStaticPools.ReturnFactorScratchDictionary(scratch); + return true; + } // private static ulong _tryFactorIntoCountsInternalHits; // private static ulong _tryCalculateCycleLengthHeuristicHits; @@ -451,471 +451,462 @@ public class MersenneDivisorCycles return true; } - private static void AddFactor(Dictionary counts, ulong factor) - { - ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists); - entry = exists ? entry + 1 : 1; - } - - private static void AddFactor(Dictionary counts, ulong factor, int multiplicity) - { - ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists); - entry = exists ? entry + multiplicity : multiplicity; - } - - private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) - { - if (factorCounts.Count == 0) - { - return initialOrder; - } - - int count = factorCounts.Count; - ulong[] primes = ThreadStaticPools.UlongPool.Rent(count); - - int index = 0; - foreach (KeyValuePair entry in factorCounts) - { - primes[index] = entry.Key; - index++; - } - - Array.Sort(primes, 0, count); - - const int StackMultiplicityCapacity = 16; - const int HeapMultiplicityCapacity = 256; - - // Scanner workloads keep almost all multiplicities within the stack-backed span, so the stack path handles the hot loop. - Span stackCandidates = stackalloc ulong[StackMultiplicityCapacity]; - Span stackEvaluations = stackalloc bool[StackMultiplicityCapacity]; - - ulong[]? heapCandidateArray = null; - bool[]? heapEvaluationArray = null; - - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - - ulong order = initialOrder; - for (int i = 0; i < count; i++) - { - ulong prime = primes[i]; - int multiplicity = factorCounts[prime]; - // Multiplicities collected by TryFactorIntoCountsInternal are strictly positive along the scanner path. - // Keeping the guard commented documents the assumption so it does not resurface in future reviews. - // if (multiplicity <= 0) - // { - // continue; - // } - - if (multiplicity <= StackMultiplicityCapacity) - { - ProcessReduceOrderPrime(stackCandidates.Slice(0, multiplicity), stackEvaluations.Slice(0, multiplicity), ref order, prime, multiplicity, ref stepper); - continue; - } - - if (multiplicity > HeapMultiplicityCapacity) - { - throw new InvalidOperationException($"Factor multiplicity {multiplicity} exceeds the supported limit of {HeapMultiplicityCapacity}."); - } - - // Larger multiplicities appear during stress tests; reuse a shared heap buffer in those cases. - heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(HeapMultiplicityCapacity); - heapEvaluationArray ??= ThreadStaticPools.BoolPool.Rent(HeapMultiplicityCapacity); - ProcessReduceOrderPrime(heapCandidateArray.AsSpan(0, multiplicity), heapEvaluationArray.AsSpan(0, multiplicity), ref order, prime, multiplicity, ref stepper); - } - - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - ThreadStaticPools.UlongPool.Return(primes); - - if (heapCandidateArray is not null) - { - ThreadStaticPools.UlongPool.Return(heapCandidateArray); - ThreadStaticPools.BoolPool.Return(heapEvaluationArray!); - } - - return order; - } - - private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper) - { - ulong working = order; - int writeIndex = multiplicity; - while (writeIndex > 0) - { - if (working < prime) - { - break; - } - - ulong reduced = working / prime; - if (reduced == 0UL) - { - break; - } - - writeIndex--; - candidateBuffer[writeIndex] = reduced; - working = reduced; - } - - int actual = multiplicity - writeIndex; - if (actual == 0) - { - return; - } - - Span candidates = candidateBuffer.Slice(writeIndex, actual); - Span evaluations = evaluationBuffer.Slice(0, actual); - - stepper.Reset(); - - evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]); - for (int j = 1; j < actual; j++) - { - evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); - } - - for (int j = actual - 1; j >= 0; j--) - { - if (!evaluations[j]) - { - break; - } - - order = candidates[j]; - } - } - - private static ulong PollardRho64(ulong n) - { - // The EvenPerfectBitScanner path removes all factors of two before invoking Pollard-Rho, so this branch never executes - // outside synthetic tests and benchmarks. Keep the guard documented to prevent future regressions without paying the hot- - // path cost. - // if ((n & 1UL) == 0UL) - // { - // return 2UL; - // } - - ulong modulus = n; - ulong nMinus1 = modulus - 1UL; - ulong nMinus2 = modulus - 2UL; - - while (true) - { - ulong c = (NextRandomUInt64() % nMinus1) + 1UL; - ulong x = (NextRandomUInt64() % nMinus2) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - x = AdvancePolynomial(x, c, modulus); - y = AdvancePolynomial(y, c, modulus); - y = AdvancePolynomial(y, c, modulus); - - ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, modulus); - } - - if (d != modulus) - { - return d; - } - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) - { - UInt128 value = (UInt128)x * x + c; - return (ulong)(value % modulus); - } - - private static ulong BinaryGcd(ulong a, ulong b) - { - if (a == 0UL) - { - return b; - } - - if (b == 0UL) - { - return a; - } - - int shift = BitOperations.TrailingZeroCount(a | b); - a >>= BitOperations.TrailingZeroCount(a); - - while (true) - { - b >>= BitOperations.TrailingZeroCount(b); - if (a > b) - { - (a, b) = (b, a); - } - - b -= a; - if (b == 0UL) - { - return a << shift; - } - } - } - - private static ulong NextRandomUInt64() - { - Span buffer = stackalloc byte[8]; - RandomNumberGenerator.Fill(buffer); - return BinaryPrimitives.ReadUInt64LittleEndian(buffer); - } - - private sealed class FactorPairValueComparer : IComparer> - { - public static readonly FactorPairValueComparer Instance = new FactorPairValueComparer(); - - private FactorPairValueComparer() - { - } - - public int Compare(KeyValuePair x, KeyValuePair y) - { - if (x.Key < y.Key) - { - return -1; - } - - if (x.Key > y.Key) - { - return 1; - } - - return 0; - } - } - - public static (long nextPosition, long completeCount) FindLast(string path) - { - using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); - long count = 0; - foreach (var _ in EnumerateStream(outputStream)) - { - count++; - } - return (outputStream.Position, count); - } - - private const int BufferSize10M = 10 * 1024 * 1024; - - public static void Generate(string path, ulong maxDivisor, int threads = 16) - { - // Binary-only generation of (divisor, cycle) pairs. CSV scaffolding removed. - ulong start = 2; - ulong blockSize = maxDivisor / (ulong)threads + 1UL; - Task[] tasks = new Task[threads]; - using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); - // Start fresh to avoid mixing formats - outputStream.SetLength(0L); - outputStream.Position = 0L; - using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); - for (var taskIndex = 0; taskIndex < threads; taskIndex++) - { - ulong threadStart = start + (ulong)taskIndex * blockSize; - if (threadStart > maxDivisor) - { - tasks[taskIndex] = Task.CompletedTask; - continue; - } - - ulong threadEnd = threadStart + blockSize - 1UL; - if (threadEnd > maxDivisor) - { - threadEnd = maxDivisor; - } - - ulong localStart = threadStart; - ulong localEnd = threadEnd; - tasks[taskIndex] = Task.Run(() => - { - ulong rangeLength = localEnd - localStart + 1UL; - (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength)); - int localCycleIndex = 0; - - try - { - for (ulong divisor = localStart; divisor <= localEnd; divisor++) - { - if ((divisor & 1UL) == 0UL) - { - continue; - } - - // TODO: Replace this modulo-based filter with the cached Mod3/Mod5/Mod7/Mod11 - // helpers so the CPU generator stops paying for `%` in the small-prime sieve - // and lines up with the divisor-cycle pipeline. - if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) - { - continue; - } - - // TODO: Migrate this generation path to the shared - // unrolled-hex calculator so the CPU generator matches - // the GPU benchmark leader for large divisors. - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - ulong cycle = CalculateCycleLength(divisor, divisorData); - localCycles[localCycleIndex++] = (divisor, cycle); - } - - if (localCycleIndex == 0) - { - return; - } - - lock (writer) - { - writer.BaseStream.Position = writer.BaseStream.Length; - for (int i = 0; i < localCycleIndex; i++) - { - (ulong divisor, ulong cycle) = localCycles[i]; - writer.Write(divisor); - writer.Write(cycle); - } - - writer.Flush(); - } - } - finally - { - ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false); - } - }); - } - - Task.WaitAll(tasks); - } - - public static void GenerateGpu(string path, ulong maxDivisor, int batchSize = 1_000_000, long skipCount = 0L, long nextPosition = 0L) - { - // Prepare output file with header - - ulong start = 3UL; - var locker = new object(); - - var pool = ArrayPool.Shared; - ulong batchSizeUL = (ulong)batchSize, d, end; - int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx; - ulong[] divisors, outCycles, validDivisors; - - divisors = pool.Rent(count); - outCycles = pool.Rent(count); - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - var accelerator = lease.Accelerator; - var stream = accelerator.CreateStream(); - var kernel = accelerator.LoadAutoGroupedStreamKernel< - Index1D, - ArrayView1D, - ArrayView1D>(DivisorCycleKernels.GpuDivisorCycleKernel); - - - using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true); - if (nextPosition > 0L) - { - outputStream.Position = nextPosition; - } - else - { - outputStream.SetLength(0L); - outputStream.Position = 0L; - } - - using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); - try - { - while (skipCount > 0L && start <= maxDivisor) - { - // TODO: Swap these modulo checks for the cached Mod3/Mod5/Mod7/Mod11 helpers so - // range initialization avoids repeated `%` calls before the unrolled-hex pipeline - // takes over. - if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL) - { - start++; - continue; - } - - start++; - skipCount--; - } - - while (start <= maxDivisor) - { - end = Math.Min(start + batchSizeUL - 1UL, maxDivisor); - count = checked((int)(end - start + 1UL)); - - idx = 0; - for (d = start; d <= end; d++) - { - if ((d & 1UL) == 0UL) - continue; - - divisors[idx++] = d; - } - - validDivisors = divisors[..idx]; - - // Use GpuKernelPool to get a kernel and context - // GpuKernelPool.Run((accelerator, stream) => - var bufferDiv = accelerator.Allocate1D(validDivisors); - var bufferCycle = accelerator.Allocate1D(idx); - - kernel( - idx, - bufferDiv.View, - bufferCycle.View); - - accelerator.Synchronize(); - - bufferCycle.View.CopyToCPU(ref outCycles[0], idx); - bufferCycle.Dispose(); - bufferDiv.Dispose(); - - // Collect results in order - for (i = 0; i < idx; i++) - { - d = divisors[i]; - - // TODO: Replace this modulo sieve with the cached Mod helpers once the divisor-cycle cache - // is mandatory so cycle enumeration can drop the slower `%` operations in this hot loop. - if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL) - continue; - - writer.Write(d); - writer.Write(outCycles[i]); - } - - writer.Flush(); - Console.WriteLine($"...processed divisor = {d}"); - start = end + 1UL; - } - } - finally - { - pool.Return(divisors, clearArray: false); - pool.Return(outCycles, clearArray: false); - } - - stream.Dispose(); - lease.Dispose(); - - } - - - public static ulong CalculateCycleLengthGpu(ulong divisor) => DivisorCycleKernels.CalculateCycleLengthGpu(divisor); - - public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false) - { - if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong cycleLength, skipPrimeOrderHeuristic)) - return cycleLength; - - return CalculateCycleLengthFallback(divisor); - } + private static void AddFactor(Dictionary counts, ulong factor) + { + ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists); + entry = exists ? entry + 1 : 1; + } + + private static void AddFactor(Dictionary counts, ulong factor, int multiplicity) + { + ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists); + entry = exists ? entry + multiplicity : multiplicity; + } + + private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) + { + if (factorCounts.Count == 0) + { + return initialOrder; + } + + int count = factorCounts.Count; + ulong[] primes = ThreadStaticPools.UlongPool.Rent(count); + + int index = 0; + foreach (KeyValuePair entry in factorCounts) + { + primes[index] = entry.Key; + index++; + } + + Array.Sort(primes, 0, count); + + const int StackMultiplicityCapacity = 16; + const int HeapMultiplicityCapacity = 256; + + // Scanner workloads keep almost all multiplicities within the stack-backed span, so the stack path handles the hot loop. + Span stackCandidates = stackalloc ulong[StackMultiplicityCapacity]; + Span stackEvaluations = stackalloc bool[StackMultiplicityCapacity]; + + ulong[]? heapCandidateArray = null; + bool[]? heapEvaluationArray = null; + + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + + ulong order = initialOrder; + for (int i = 0; i < count; i++) + { + ulong prime = primes[i]; + int multiplicity = factorCounts[prime]; + // Multiplicities collected by TryFactorIntoCountsInternal are strictly positive along the scanner path. + // Keeping the guard commented documents the assumption so it does not resurface in future reviews. + // if (multiplicity <= 0) + // { + // continue; + // } + + if (multiplicity <= StackMultiplicityCapacity) + { + ProcessReduceOrderPrime(stackCandidates.Slice(0, multiplicity), stackEvaluations.Slice(0, multiplicity), ref order, prime, multiplicity, ref stepper); + continue; + } + + if (multiplicity > HeapMultiplicityCapacity) + { + throw new InvalidOperationException($"Factor multiplicity {multiplicity} exceeds the supported limit of {HeapMultiplicityCapacity}."); + } + + // Larger multiplicities appear during stress tests; reuse a shared heap buffer in those cases. + heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(HeapMultiplicityCapacity); + heapEvaluationArray ??= ThreadStaticPools.BoolPool.Rent(HeapMultiplicityCapacity); + ProcessReduceOrderPrime(heapCandidateArray.AsSpan(0, multiplicity), heapEvaluationArray.AsSpan(0, multiplicity), ref order, prime, multiplicity, ref stepper); + } + + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + ThreadStaticPools.UlongPool.Return(primes); + + if (heapCandidateArray is not null) + { + ThreadStaticPools.UlongPool.Return(heapCandidateArray); + ThreadStaticPools.BoolPool.Return(heapEvaluationArray!); + } + + return order; + } + + private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper) + { + ulong working = order; + int writeIndex = multiplicity; + while (writeIndex > 0) + { + if (working < prime) + { + break; + } + + ulong reduced = working / prime; + if (reduced == 0UL) + { + break; + } + + writeIndex--; + candidateBuffer[writeIndex] = reduced; + working = reduced; + } + + int actual = multiplicity - writeIndex; + if (actual == 0) + { + return; + } + + Span candidates = candidateBuffer.Slice(writeIndex, actual); + Span evaluations = evaluationBuffer.Slice(0, actual); + + stepper.Reset(); + + evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]); + for (int j = 1; j < actual; j++) + { + evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); + } + + for (int j = actual - 1; j >= 0; j--) + { + if (!evaluations[j]) + { + break; + } + + order = candidates[j]; + } + } + + private static ulong PollardRho64(ulong n) + { + // The EvenPerfectBitScanner path removes all factors of two before invoking Pollard-Rho, so this branch never executes + // outside synthetic tests and benchmarks. Keep the guard documented to prevent future regressions without paying the hot- + // path cost. + // if ((n & 1UL) == 0UL) + // { + // return 2UL; + // } + + ulong modulus = n; + ulong nMinus1 = modulus - 1UL; + ulong nMinus2 = modulus - 2UL; + + while (true) + { + ulong c = (NextRandomUInt64() % nMinus1) + 1UL; + ulong x = (NextRandomUInt64() % nMinus2) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomial(x, c, modulus); + y = AdvancePolynomial(y, c, modulus); + y = AdvancePolynomial(y, c, modulus); + + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, modulus); + } + + if (d != modulus) + { + return d; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) + { + UInt128 value = (UInt128)x * x + c; + return (ulong)(value % modulus); + } + + private static ulong BinaryGcd(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + a >>= BitOperations.TrailingZeroCount(a); + + while (true) + { + b >>= BitOperations.TrailingZeroCount(b); + if (a > b) + { + (a, b) = (b, a); + } + + b -= a; + if (b == 0UL) + { + return a << shift; + } + } + } + + private static ulong NextRandomUInt64() + { + Span buffer = stackalloc byte[8]; + RandomNumberGenerator.Fill(buffer); + return BinaryPrimitives.ReadUInt64LittleEndian(buffer); + } + + private sealed class FactorPairValueComparer : IComparer> + { + public static readonly FactorPairValueComparer Instance = new FactorPairValueComparer(); + + private FactorPairValueComparer() + { + } + + public int Compare(KeyValuePair x, KeyValuePair y) + { + if (x.Key < y.Key) + { + return -1; + } + + if (x.Key > y.Key) + { + return 1; + } + + return 0; + } + } + + public static (long nextPosition, long completeCount) FindLast(string path) + { + using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); + long count = 0; + foreach (var _ in EnumerateStream(outputStream)) + { + count++; + } + return (outputStream.Position, count); + } + + private const int BufferSize10M = 10 * 1024 * 1024; + + public static void Generate(string path, ulong maxDivisor, int threads = 16) + { + // Binary-only generation of (divisor, cycle) pairs. CSV scaffolding removed. + ulong start = 2; + ulong blockSize = maxDivisor / (ulong)threads + 1UL; + Task[] tasks = new Task[threads]; + using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); + // Start fresh to avoid mixing formats + outputStream.SetLength(0L); + outputStream.Position = 0L; + using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); + for (var taskIndex = 0; taskIndex < threads; taskIndex++) + { + ulong threadStart = start + (ulong)taskIndex * blockSize; + if (threadStart > maxDivisor) + { + tasks[taskIndex] = Task.CompletedTask; + continue; + } + + ulong threadEnd = threadStart + blockSize - 1UL; + if (threadEnd > maxDivisor) + { + threadEnd = maxDivisor; + } + + ulong localStart = threadStart; + ulong localEnd = threadEnd; + tasks[taskIndex] = Task.Run(() => + { + ulong rangeLength = localEnd - localStart + 1UL; + (ulong divisor, ulong cycle)[] localCycles = ThreadStaticPools.RentUlongUlongArray(checked((int)rangeLength)); + // (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength)); + int localCycleIndex = 0; + + for (ulong divisor = localStart; divisor <= localEnd; divisor++) + { + if ((divisor & 1UL) == 0UL) + { + continue; + } + + // TODO: Replace this modulo-based filter with the cached Mod3/Mod5/Mod7/Mod11 + // helpers so the CPU generator stops paying for `%` in the small-prime sieve + // and lines up with the divisor-cycle pipeline. + if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) + { + continue; + } + + // TODO: Migrate this generation path to the shared + // unrolled-hex calculator so the CPU generator matches + // the GPU benchmark leader for large divisors. + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong cycle = CalculateCycleLength(divisor, divisorData); + localCycles[localCycleIndex++] = (divisor, cycle); + } + + if (localCycleIndex == 0) + { + ThreadStaticPools.Return(localCycles); + return; + } + + lock (writer) + { + writer.BaseStream.Position = writer.BaseStream.Length; + for (int i = 0; i < localCycleIndex; i++) + { + (ulong divisor, ulong cycle) = localCycles[i]; + writer.Write(divisor); + writer.Write(cycle); + } + + writer.Flush(); + } + + ThreadStaticPools.Return(localCycles); + }); + } + + Task.WaitAll(tasks); + } + + public static void GenerateGpu(string path, ulong maxDivisor, int batchSize = 1_000_000, long skipCount = 0L, long nextPosition = 0L) + { + // Prepare output file with header + + ulong start = 3UL; + var locker = new object(); + + var pool = ArrayPool.Shared; + ulong batchSizeUL = (ulong)batchSize, d, end; + int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx; + ulong[] divisors, outCycles, validDivisors; + + divisors = pool.Rent(count); + outCycles = pool.Rent(count); + var lease = GpuKernelPool.GetKernel(); + var accelerator = lease.Accelerator; + var stream = lease.Stream; + var kernel = accelerator.LoadAutoGroupedStreamKernel< + Index1D, + ArrayView1D, + ArrayView1D>(DivisorCycleKernels.GpuDivisorCycleKernel); + + + using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true); + if (nextPosition > 0L) + { + outputStream.Position = nextPosition; + } + else + { + outputStream.SetLength(0L); + outputStream.Position = 0L; + } + + using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); + while (skipCount > 0L && start <= maxDivisor) + { + // TODO: Swap these modulo checks for the cached Mod3/Mod5/Mod7/Mod11 helpers so + // range initialization avoids repeated `%` calls before the unrolled-hex pipeline + // takes over. + if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL) + { + start++; + continue; + } + + start++; + skipCount--; + } + + while (start <= maxDivisor) + { + end = Math.Min(start + batchSizeUL - 1UL, maxDivisor); + count = checked((int)(end - start + 1UL)); + + idx = 0; + for (d = start; d <= end; d++) + { + if ((d & 1UL) == 0UL) + continue; + + divisors[idx++] = d; + } + + validDivisors = divisors[..idx]; + + // Use GpuKernelPool to get a kernel and context + // GpuKernelPool.Run((accelerator, stream) => + var bufferDiv = accelerator.Allocate1D(validDivisors); + var bufferCycle = accelerator.Allocate1D(idx); + + kernel( + idx, + bufferDiv.View, + bufferCycle.View); + + stream.Synchronize(); + + bufferCycle.View.CopyToCPU(stream, ref outCycles[0], idx); + bufferCycle.Dispose(); + bufferDiv.Dispose(); + + // Collect results in order + for (i = 0; i < idx; i++) + { + d = divisors[i]; + + // TODO: Replace this modulo sieve with the cached Mod helpers once the divisor-cycle cache + // is mandatory so cycle enumeration can drop the slower `%` operations in this hot loop. + if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL) + continue; + + writer.Write(d); + writer.Write(outCycles[i]); + } + + writer.Flush(); + Console.WriteLine($"...processed divisor = {d}"); + start = end + 1UL; + } + + pool.Return(divisors, clearArray: false); + pool.Return(outCycles, clearArray: false); + + stream.Dispose(); + lease.Dispose(); + } + + + public static ulong CalculateCycleLengthGpu(ulong divisor) => DivisorCycleKernels.CalculateCycleLengthGpu(divisor); + + public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false) + { + if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong cycleLength, skipPrimeOrderHeuristic)) + return cycleLength; + + return CalculateCycleLengthFallback(divisor); + } // [ThreadStatic] // private static HeuristicPrimeTester? _primeTester; @@ -962,18 +953,18 @@ public class MersenneDivisorCycles return false; } - private static ulong CalculateCycleLengthFallback(ulong divisor) - { - // Otherwise, find order of 2 mod divisor using the shared unrolled-hex stepping sequence. - return divisor.CalculateMersenneDivisorCycleLengthUnrolledHexCpu(); - } + private static ulong CalculateCycleLengthFallback(ulong divisor) + { + // Otherwise, find order of 2 mod divisor using the shared unrolled-hex stepping sequence. + return divisor.CalculateMersenneDivisorCycleLengthUnrolledHexCpu(); + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void EnsureSmallBuffer() - { - if (_smallCycles is null) - { - _smallCycles = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; - } - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void EnsureSmallBuffer() + { + if (_smallCycles is null) + { + _smallCycles = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + } + } } diff --git a/PerfectNumbers.Core/MersenneNumberTester.cs b/PerfectNumbers.Core/MersenneNumberTester.cs index de77340..a8620a3 100644 --- a/PerfectNumbers.Core/MersenneNumberTester.cs +++ b/PerfectNumbers.Core/MersenneNumberTester.cs @@ -146,7 +146,7 @@ public sealed class MersenneNumberTester( UInt128[] qs = qsBuffer; - var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); + var gpuLease = GpuKernelPool.GetKernel(); var accelerator = gpuLease.Accelerator; var stream = gpuLease.Stream; var orderKernel = gpuLease.OrderKernel!; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 126603d..fa3120d 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,8 +1,6 @@ -using System; using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; -using ILGPU.Runtime.OpenCL; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -21,40 +19,53 @@ internal static partial class PrimeOrderCalculator Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); remaining = value; - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var lease = GpuKernelPool.GetKernel(); Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; + var stream = lease.Stream; + KernelContainer kernels = lease.Kernels; - SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(lease.Kernels, accelerator); - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(lease.Kernels, accelerator, GpuSmallPrimeFactorSlots); + ScratchBuffer scratch = GpuScratchBufferPool.Rent(accelerator, GpuSmallPrimeFactorSlots, 0); + + SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(kernels, accelerator); + + if (kernels.SmallPrimeFactorsPrimes is null) throw new NullReferenceException($"{nameof(kernels.SmallPrimeFactorsPrimes)} is null"); + if (kernels.SmallPrimeFactorsSquares is null) throw new NullReferenceException($"{nameof(kernels.SmallPrimeFactorsSquares)} is null"); var kernel = lease.SmallPrimeFactorKernel; + + ArrayView1D smallPrimeFactorCountSlotView = scratch.SmallPrimeFactorCountSlot.View; + ArrayView1D smallPrimeFactorRemainingSlotView = scratch.SmallPrimeFactorRemainingSlot.View; + ArrayView1D smallPrimeFactorPrimeSlotsView = scratch.SmallPrimeFactorPrimeSlots.View; + ArrayView1D smallPrimeFactorExponentSlotsView = scratch.SmallPrimeFactorExponentSlots.View; + ArrayView1D smallPrimeFactorsPrimesView = kernels.SmallPrimeFactorsPrimes!.View; + kernel( stream, 1, value, limit, - tables.PrimesView, - tables.SquaresView, - tables.Count, - scratch.PrimeSlots.View, - scratch.ExponentSlots.View, - scratch.CountSlot.View, - scratch.RemainingSlot.View); + smallPrimeFactorsPrimesView, + kernels.SmallPrimeFactorsSquares!.View, + (int)smallPrimeFactorsPrimesView.Length, + smallPrimeFactorPrimeSlotsView, + smallPrimeFactorExponentSlotsView, + smallPrimeFactorCountSlotView, + smallPrimeFactorRemainingSlotView); stream.Synchronize(); factorCount = 0; - scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); + smallPrimeFactorCountSlotView.CopyToCPU(stream, 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); + smallPrimeFactorPrimeSlotsView.CopyToCPU(stream, ref MemoryMarshal.GetReference(primeBuffer), factorCount); + + smallPrimeFactorExponentSlotsView.CopyToCPU(stream, ref MemoryMarshal.GetReference(exponentBuffer), factorCount); } - scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); + smallPrimeFactorRemainingSlotView.CopyToCPU(stream, ref remaining, 1); for (int i = 0; i < factorCount; i++) { @@ -63,9 +74,8 @@ internal static partial class PrimeOrderCalculator counts.Add(primeValue, exponent); } + GpuScratchBufferPool.Return(scratch); lease.Dispose(); - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); return true; } @@ -88,32 +98,40 @@ internal static partial class PrimeOrderCalculator factorSpan[i] = factors[i].Value; } - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + GpuKernelLease lease = GpuKernelPool.GetKernel(); Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - SpecialMaxScratch scratch = GpuKernelPool.EnsureSpecialMaxScratch(lease.Kernels, accelerator, factorCapacity: factorCount); + var stream = lease.Stream; - scratch.FactorsView.SubView(0, factorCount).CopyFromCPU(ref MemoryMarshal.GetReference(factorSpan), factorCount); + ScratchBuffer scratch = GpuScratchBufferPool.Rent(accelerator, 0, factorCount); + + ArrayView1D specialMaxFactorsView = scratch.SpecialMaxFactors.View; + specialMaxFactorsView.SubView(0, factorCount).CopyFromCPU(stream, ref MemoryMarshal.GetReference(factorSpan), factorCount); var kernel = lease.SpecialMaxKernel; + + ArrayView1D specialMaxResultView = scratch.SpecialMaxResult.View; + kernel( stream, 1, phi, - scratch.FactorsView, + specialMaxFactorsView, factorCount, divisorData, - scratch.CandidatesView, - scratch.ResultView); - + scratch.SpecialMaxCandidates.View, + specialMaxResultView); stream.Synchronize(); - ushort result = 0; - scratch.ResultView.CopyToCPU(ref result, 1); - lease.Dispose(); + Span result = stackalloc ulong[1]; + // specialMaxResultView.CopyToCPU(stream, ref result, 1); + specialMaxResultView.CopyToCPU(stream, in result); - return result != 0; + // stream.Synchronize(); + // return result != 0; + GpuScratchBufferPool.Return(scratch); + lease.Dispose(); + return result[0] != 0; } } diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 6e64b49..96e02f3 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -98,6 +98,7 @@ public sealed class PrimeTester GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTesterGpuContextPool.Rent(1); var accelerator = gpu.Accelerator; + var stream = gpu.Stream; var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.Kernel; ulong value = n; @@ -108,7 +109,7 @@ public sealed class PrimeTester var inputView = input.View; var outputView = output.View; - inputView.CopyFromCPU(ref value, 1); + inputView.CopyFromCPU(stream, ref value, 1); kernel( 1, inputView, @@ -123,8 +124,8 @@ public sealed class PrimeTester gpu.DevicePrimesPow2LastThree.View, gpu.DevicePrimesPow2LastNine.View, outputView); - accelerator.Synchronize(); - outputView.CopyToCPU(ref flag, 1); + stream.Synchronize(); + outputView.CopyToCPU(stream, ref flag, 1); gpu.Dispose(); GpuPrimeWorkLimiter.Release(); @@ -163,6 +164,7 @@ public sealed class PrimeTester GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTesterGpuContextPool.Rent(GpuBatchSize); var accelerator = gpu.Accelerator; + var stream = gpu.Stream; var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.Kernel; int totalLength = values.Length; int batchSize = GpuBatchSize; @@ -196,10 +198,11 @@ public sealed class PrimeTester gpu.DevicePrimesPow2LastThree.View, gpu.DevicePrimesPow2LastNine.View, outputView); - accelerator.Synchronize(); + stream.Synchronize(); + var resultSlice = results.Slice(pos, count); ref byte resultRef = ref MemoryMarshal.GetReference(resultSlice); - outputView.CopyToCPU(ref resultRef, count); + outputView.CopyToCPU(stream, ref resultRef, count); pos += count; } @@ -266,8 +269,9 @@ public sealed class PrimeTester return; } - KernelContainer kernels = GpuKernelPool.GetKernels(); - GpuStaticTableInitializer.EnsureStaticTables(kernels, SharedGpuContext.Accelerator); + Accelerator accelerator = SharedGpuContext.Accelerator; + KernelContainer kernels = GpuKernelPool.GetKernels(accelerator); + GpuStaticTableInitializer.EnsureStaticTables(kernels, accelerator); int toCreate = target - WarmedLeaseCount; for (int i = 0; i < toCreate; i++) @@ -275,7 +279,6 @@ public sealed class PrimeTester Pool.Enqueue(new PrimeTesterGpuContextLease(minBufferCapacity)); } - kernels.Dispose(); WarmedLeaseCount = target; } } @@ -379,6 +382,21 @@ public sealed class PrimeTester internal readonly Context AcceleratorContext; public readonly Accelerator Accelerator; + private AcceleratorStream? _stream; + public AcceleratorStream Stream + { + get + { + if (_stream is { } stream) + { + return stream; + } + + stream = Accelerator.CreateStream(); + _stream = stream; + return stream; + } + } public static readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> Kernel = SharedGpuContext.Accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); @@ -445,6 +463,7 @@ public sealed class PrimeTester Input?.Dispose(); Output?.Dispose(); HeuristicFlag?.Dispose(); + Input = Accelerator.Allocate1D(minCapacity); Output = Accelerator.Allocate1D(minCapacity); HeuristicFlag = Accelerator.Allocate1D(minCapacity); @@ -464,6 +483,12 @@ public sealed class PrimeTester } _returned = true; + if (_stream is { } stream) + { + stream.Dispose(); + _stream = null; + } + Return(this); } @@ -734,6 +759,7 @@ public sealed class PrimeTester // and divisor-cycle staging instead of allocating new device buffers per call. var gpu = PrimeTesterGpuContextPool.Rent(); var accelerator = gpu.Accelerator; + var stream = gpu.Stream; var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel); int length = values.Length; @@ -743,10 +769,10 @@ public sealed class PrimeTester ArrayPool pool = ThreadStaticPools.UlongPool; ulong[] temp = pool.Rent(length); values.CopyTo(temp); - inputBuffer.View.CopyFromCPU(ref temp[0], length); + inputBuffer.View.CopyFromCPU(stream, ref temp[0], length); kernel(length, inputBuffer.View, resultBuffer.View); - accelerator.Synchronize(); - resultBuffer.View.CopyToCPU(ref results[0], length); + stream.Synchronize(); + resultBuffer.View.CopyToCPU(stream, ref results[0], length); pool.Return(temp, clearArray: false); resultBuffer.Dispose(); inputBuffer.Dispose(); diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 3d5ca9b..29a7bbe 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -180,6 +180,20 @@ namespace PerfectNumbers.Core } [ThreadStatic] + private static ArrayPool<(ulong, ulong)>? _ulongUlongTupleArrayPool; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static (ulong, ulong)[] RentUlongUlongArray(int capacityHint) + { + _ulongUlongTupleArrayPool ??= ArrayPool<(ulong, ulong)>.Create(); + ArrayPool<(ulong, ulong)>? pool = _ulongUlongTupleArrayPool; + return pool.Rent(capacityHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void Return((ulong, ulong)[] array) => _ulongUlongTupleArrayPool!.Return(array, clearArray: false); + + [ThreadStatic] private static List>? _primeOrderPendingEntryListPool; [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index d4fc788..85993ae 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -456,7 +456,7 @@ public static partial class ULongExtensions { ulong result = 0UL; - GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true); + GpuKernelLease lease = GpuKernelPool.GetKernel(); Accelerator accelerator = lease.Accelerator; // Keep this commented out. It should never happen in production code. @@ -479,7 +479,7 @@ public static partial class ULongExtensions kernel(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); stream.Synchronize(); - resultBuffer.View.CopyToCPU(ref result, 1); + resultBuffer.View.CopyToCPU(stream, ref result, 1); resultBuffer.Dispose(); exponentBuffer.Dispose();