diff --git a/PerfectNumbers.Core/DeterministicRandom.cs b/PerfectNumbers.Core/Cpu/DeterministicRandomCpu.cs similarity index 92% rename from PerfectNumbers.Core/DeterministicRandom.cs rename to PerfectNumbers.Core/Cpu/DeterministicRandomCpu.cs index 279da51..3d04d1b 100644 --- a/PerfectNumbers.Core/DeterministicRandom.cs +++ b/PerfectNumbers.Core/Cpu/DeterministicRandomCpu.cs @@ -1,6 +1,6 @@ -namespace PerfectNumbers.Core; +namespace PerfectNumbers.Core.Cpu; -internal static class DeterministicRandom +internal static class DeterministicRandomCpu { [ThreadStatic] private static ulong s_threadState; diff --git a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs new file mode 100644 index 0000000..f1cb530 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs @@ -0,0 +1,36 @@ +using System.Collections.Concurrent; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu.Accelerators; + +public static class AcceleratorStreamPool +{ + private static readonly ConcurrentDictionary> _streams = new(20_480, PerfectNumberConstants.RollingAccelerators); + // private static readonly ConcurrentDictionary _streams = new(AcceleratorPool.Shared.Accelerators.ToDictionary(static x => x, CreateStream)); + + private static SemaphoreSlim CreateLock() => new(PerfectNumberConstants.ThreadsByStream); + private static AcceleratorStream CreateStream(Accelerator accelerator) => accelerator.CreateStream(); + + private static readonly ConcurrentDictionary _locks = new(AcceleratorPool.Shared.Accelerators.ToDictionary(x => x, static x => CreateLock())); + + public static AcceleratorStream Rent(Accelerator accelerator) + { + var streamLock = _locks.GetOrAdd(accelerator, static _ => new SemaphoreSlim(PerfectNumberConstants.ThreadsByStream)); + + streamLock.Wait(); + var queue = _streams.GetOrAdd(accelerator, static _ => new()); + + return queue.TryDequeue(out var stream) + ? stream + : accelerator.CreateStream(); + } + + public static void Return(AcceleratorStream stream) + { + Accelerator accelerator = stream.Accelerator; + var streamLock = _locks.GetOrAdd(accelerator, static _ => CreateLock()); + var queue = _streams.GetOrAdd(accelerator, static _ => new()); + queue.Enqueue(stream); + streamLock.Release(); + } +} diff --git a/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs index ce25c83..13a913a 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs @@ -12,7 +12,7 @@ public sealed class Pow2MontgomeryAccelerator [ThreadStatic] private static Queue? _pool; - private static AcceleratorPool _accelerators = new(PerfectNumberConstants.RollingAccelerators << 1); + private static AcceleratorPool _accelerators = new(PerfectNumberConstants.RollingAccelerators); [MethodImpl(MethodImplOptions.AggressiveInlining)] public static Pow2MontgomeryAccelerator Rent(int primeTesterCapacity) @@ -25,8 +25,8 @@ public sealed class Pow2MontgomeryAccelerator } else { - gpu.Stream = gpu.Accelerator.CreateStream(); - gpu.EnsureCapacity(PerfectNumberConstants.DefaultFactorsBuffer,primeTesterCapacity); + // gpu.Stream = gpu.Accelerator.CreateStream(); + gpu.EnsureCapacity(PerfectNumberConstants.DefaultFactorsBuffer, primeTesterCapacity); } return gpu; @@ -35,8 +35,8 @@ public sealed class Pow2MontgomeryAccelerator [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void Return(Pow2MontgomeryAccelerator gpu) { - gpu.Stream!.Dispose(); - gpu.Stream = null; + // gpu.Stream!.Dispose(); + // gpu.Stream = null; _pool!.Enqueue(gpu); } @@ -87,9 +87,9 @@ public sealed class Pow2MontgomeryAccelerator public readonly Accelerator Accelerator; public readonly MemoryBuffer1D Input; - public readonly MemoryBuffer1D Output; + public readonly MemoryBuffer1D OutputUlong; public MemoryBuffer1D PrimeTestInput; - public MemoryBuffer1D PrimeTestOutput; + public MemoryBuffer1D OutputByte; public MemoryBuffer1D HeuristicFlag; public MemoryBuffer1D, Stride1D.Dense> Pow2ModEntriesToTestOnDevice; public Dictionary Pow2ModEntriesToTestOnHost = []; @@ -101,14 +101,14 @@ public sealed class Pow2MontgomeryAccelerator public readonly MemoryBuffer1D DevicePrimesPow2LastSeven; public readonly MemoryBuffer1D DevicePrimesPow2LastThree; public readonly MemoryBuffer1D DevicePrimesPow2LastNine; - public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> SmallPrimeSieveKernel; + public readonly Kernel SmallPrimeSieveKernel; - public AcceleratorStream? Stream; - - public readonly Action, MontgomeryDivisorData, ArrayView1D> ConvertToStandardKernel; - public readonly Action, MontgomeryDivisorData, ArrayView1D> KeepMontgomeryKernel; + // public AcceleratorStream? Stream; public readonly Kernel CheckFactorsKernel; + public readonly Kernel ConvertToStandardKernel; + public readonly Kernel KeepMontgomeryKernel; + public readonly Kernel PollardRhoKernel; [ThreadStatic] private static Dictionary, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>? SmallPrimeSieveKernelPool; @@ -137,15 +137,14 @@ public sealed class Pow2MontgomeryAccelerator if (primeTesterCapacity > PrimeTestInput.Length) { - PrimeTestInput?.Dispose(); - PrimeTestOutput?.Dispose(); + OutputByte?.Dispose(); HeuristicFlag?.Dispose(); // lock (Accelerator) { PrimeTestInput = Accelerator.Allocate1D(primeTesterCapacity); - PrimeTestOutput = Accelerator.Allocate1D(primeTesterCapacity); + OutputByte = Accelerator.Allocate1D(primeTesterCapacity); HeuristicFlag = Accelerator.Allocate1D(primeTesterCapacity); } } @@ -154,14 +153,16 @@ public sealed class Pow2MontgomeryAccelerator public Pow2MontgomeryAccelerator(Accelerator accelerator, int factorsCount, int primeTesterCapacity) { Accelerator = accelerator; - AcceleratorStream stream = accelerator.CreateStream(); - Stream = stream; + AcceleratorStream stream = AcceleratorStreamPool.Rent(accelerator); + // AcceleratorStream stream = accelerator.DefaultStream; + // AcceleratorStream stream = accelerator.CreateStream(); + // Stream = stream; Input = accelerator.Allocate1D(1); - Output = accelerator.Allocate1D(1); + OutputUlong = accelerator.Allocate1D(1); Pow2ModEntriesToTestOnDevice = accelerator.Allocate1D>(factorsCount); PrimeTestInput = Accelerator.Allocate1D(primeTesterCapacity); - PrimeTestOutput = Accelerator.Allocate1D(primeTesterCapacity); + OutputByte = Accelerator.Allocate1D(primeTesterCapacity); HeuristicFlag = Accelerator.Allocate1D(primeTesterCapacity); var sharedTables = LastDigitGpuTables.EnsureStaticTables(accelerator, stream); @@ -174,30 +175,30 @@ public sealed class Pow2MontgomeryAccelerator DevicePrimesPow2LastThree = sharedTables.DevicePrimesPow2LastThree; DevicePrimesPow2LastNine = sharedTables.DevicePrimesPow2LastNine; - var keepKernel = accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery); - var keepLauncher = KernelUtil.GetKernel(keepKernel); - - KeepMontgomeryKernel = keepLauncher.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); - - var convertKernel = accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelConvertToStandard); - var convertLauncher = KernelUtil.GetKernel(convertKernel); - - ConvertToStandardKernel = convertLauncher.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); - var checkFactorsKernel = accelerator.LoadStreamKernel, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D>(PrimeOrderGpuHeuristics.CheckFactorsKernel); CheckFactorsKernel = KernelUtil.GetKernel(checkFactorsKernel); - SmallPrimeSieveKernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel)).CreateLauncherDelegate, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>(); + var convertKernel = accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelConvertToStandard); + ConvertToStandardKernel = KernelUtil.GetKernel(convertKernel); + + var keepKernel = accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery); + KeepMontgomeryKernel = KernelUtil.GetKernel(keepKernel); + + SmallPrimeSieveKernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel)); + + PollardRhoKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, ArrayView1D, ArrayView1D>(Pow2MontgomeryKernels.TryPollardRhoKernel)); + + AcceleratorStreamPool.Return(stream); } public void Dispose() { Input.Dispose(); - Output.Dispose(); - Stream?.Dispose(); + OutputUlong.Dispose(); + // Stream?.Dispose(); Pow2ModEntriesToTestOnDevice.Dispose(); PrimeTestInput?.Dispose(); - PrimeTestOutput?.Dispose(); + OutputByte?.Dispose(); HeuristicFlag?.Dispose(); // Stream is disposed upon return to the queue, so we don't need to dispose it here. diff --git a/PerfectNumbers.Core/Gpu/DeterministicRandomGpu.cs b/PerfectNumbers.Core/Gpu/DeterministicRandomGpu.cs new file mode 100644 index 0000000..996ba67 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/DeterministicRandomGpu.cs @@ -0,0 +1,75 @@ +using System.Runtime.CompilerServices; +using ILGPU; + +namespace PerfectNumbers.Core.Gpu; + +public static class ThreadStaticDeterministicRandomGpu +{ + [ThreadStatic] + private static DeterministicRandomGpu? _exclusive; + internal static DeterministicRandomGpu Exclusive => _exclusive ?? new DeterministicRandomGpu(); +} + +internal struct DeterministicRandomGpu +{ + private static long s_seedCounter; + + private ulong _state; + public readonly ulong State + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _state; + } + } + + public DeterministicRandomGpu() + { + _state = CreateThreadState(); + } + + public DeterministicRandomGpu(ulong initialState) + { + _state = initialState; + } + + public ulong NextUInt64() + { + ulong state = _state; + state ^= state >> 12; + state ^= state << 25; + state ^= state >> 27; + _state = state; + return state * 2685821657736338717UL; + } + + public GpuUInt128 NextUInt128() + { + GpuUInt128 high = (GpuUInt128)NextUInt64(); + ReadOnlyGpuUInt128 low = (ReadOnlyGpuUInt128)NextUInt64(); + high.ShiftLeft(64); + high.Add(low); + return high; + } + + public void SetState(ulong state) + { + _state = state; + } + + private static ulong CreateThreadState() + { + ulong seed = (ulong)Atomic.Add(ref s_seedCounter, 1); + seed += 0x9E3779B97F4A7C15UL; + seed = (seed ^ (seed >> 30)) * 0xBF58476D1CE4E5B9UL; + seed = (seed ^ (seed >> 27)) * 0x94D049BB133111EBUL; + seed ^= seed >> 31; + if (seed == 0UL) + { + seed = 0x9E3779B97F4A7C15UL; + } + + return seed; + } +} diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index 5159783..b134fad 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -261,13 +261,13 @@ public struct GpuUInt128 : IComparable, IEquatable public override readonly int GetHashCode() => HashCode.Combine(High, Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Add(in ReadOnlyGpuUInt128 other) => AddInternal(other.High, other.Low); + public void Add(in ReadOnlyGpuUInt128 other) => Add(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Add(GpuUInt128 other) => AddInternal(other.High, other.Low); + public void Add(GpuUInt128 other) => Add(other.High, other.Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void AddInternal(ulong otherHigh, ulong otherLow) + public void Add(ulong otherHigh, ulong otherLow) { ulong originalLow = Low; // We're reusing originalLow for carry to avoid an extra local. @@ -316,7 +316,7 @@ public struct GpuUInt128 : IComparable, IEquatable return; } - AddInternal(valueHigh, valueLow); + Add(valueHigh, valueLow); if (CompareTo(in modulus) >= 0) { Sub(in modulus); @@ -589,9 +589,6 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 BinaryGcd(GpuUInt128 u, GpuUInt128 v) { - // TODO: Replace this scalar binary GCD with the branchless reduction from - // GpuUInt128BinaryGcdBenchmarks so CPU fallbacks stay aligned with the GPU kernel - // performance when resolving large divisor residues. if (u.IsZero) { return v; @@ -602,26 +599,50 @@ public struct GpuUInt128 : IComparable, IEquatable return u; } - int shift = TrailingZeroCount(new GpuUInt128(u.High | v.High, u.Low | v.Low)); + ulong combinedLow = u.Low | v.Low; + int shift; + if (combinedLow != 0UL) + { + shift = BitOperations.TrailingZeroCount(combinedLow); + } + else + { + ulong combinedHigh = u.High | v.High; + shift = 64 + BitOperations.TrailingZeroCount(combinedHigh); + } + int zu = TrailingZeroCount(u); u >>= zu; - do + GpuUInt128 currentV = v; + while (true) { - int zv = TrailingZeroCount(v); - v >>= zv; - if (u > v) + int zv = TrailingZeroCount(currentV); + currentV >>= zv; + + bool swap = u > currentV; + GpuUInt128 min = ConditionalSelect(u, currentV, swap); + GpuUInt128 max = ConditionalSelect(currentV, u, swap); + max.Sub(min); + + if (max.IsZero) { - (u, v) = (v, u); + return min << shift; } - v -= u; + u = min; + currentV = max; } - while (!v.IsZero); - - return u << shift; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static GpuUInt128 ConditionalSelect(in GpuUInt128 left, in GpuUInt128 right, bool useRight) + { + ulong mask = useRight ? ulong.MaxValue : 0UL; + return new GpuUInt128( + left.High ^ ((left.High ^ right.High) & mask), + left.Low ^ ((left.Low ^ right.Low) & mask)); + } public void Sub(ulong value) { @@ -1529,6 +1550,15 @@ public struct GpuUInt128 : IComparable, IEquatable high = hi; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Pow2() + { + ulong operand = Low; + Low = operand * operand; + ulong highProduct = operand * High; + High = highProduct + MulHigh(operand, operand); + } + } // TODO: Check if the TODO below is still relevant. diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs index 306b2b8..3ac3631 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs @@ -1,3 +1,4 @@ +using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; @@ -5,21 +6,75 @@ namespace PerfectNumbers.Core.Gpu.Kernels; internal static class Pow2MontgomeryKernels { - public static void Pow2MontgomeryKernelKeepMontgomery( - Index1D index, - ArrayView1D exponents, - MontgomeryDivisorData divisor, - ArrayView1D results) - { - results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponents[index]); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static GpuUInt128 AdvancePolynomialGpu(GpuUInt128 x, ulong c, ulong modulus) + { + x.Pow2(); + x.AddMod(c, modulus); + return x; + } - public static void Pow2MontgomeryKernelConvertToStandard( - Index1D index, - ArrayView1D exponents, - MontgomeryDivisorData divisor, - ArrayView1D results) - { - results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponents[index]); - } + public static void Pow2MontgomeryKernelKeepMontgomery( + Index1D index, + ArrayView1D exponents, + MontgomeryDivisorData divisor, + ArrayView1D results) + { + results[0] = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponents[0]); + } + + public static void Pow2MontgomeryKernelConvertToStandard( + Index1D index, + ArrayView1D exponents, + MontgomeryDivisorData divisor, + ArrayView1D results) + { + results[0] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponents[0]); + } + + + public static void TryPollardRhoKernel(ulong n, int limit, ArrayView1D randomState, ArrayView1D factored, ArrayView1D factor) + { + factored[0] = 0; + if ((n & 1UL) == 0UL) + { + factor[0] = 2UL; + factored[0] = 1; + return; + } + + int attempt = 0; + GpuUInt128 nGpuUInt128 = (GpuUInt128)n; + DeterministicRandomGpu random = new(randomState[0]); + while (attempt < limit) + { + ulong c = (random.NextUInt64() % (n - 1UL)) + 1UL; + GpuUInt128 x = (random.NextUInt64() % (n - 2UL)) + 2UL; + GpuUInt128 y = x; + GpuUInt128 d = GpuUInt128.One; + + while (d.CompareTo(GpuUInt128.One) == 0) + { + x = AdvancePolynomialGpu(x, c, modulus: n); + y = AdvancePolynomialGpu(y, c, n); + y = AdvancePolynomialGpu(y, c, n); + GpuUInt128 diff = x.CompareTo(y) > 0 ? x - y : y - x; + d = GpuUInt128.BinaryGcd(diff, nGpuUInt128); + } + + if (d != n) + { + if (d.High == 0) + { + factor[0] = d.Low; + factored[0] = 1; + return; + } + } + + attempt++; + } + + randomState[0] = random.State; + } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs index 3ba4077..17ab513 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs @@ -2,6 +2,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core.Cpu; namespace PerfectNumbers.Core.Gpu; @@ -1410,5 +1411,4 @@ internal static partial class PrimeOrderGpuHeuristics ulong directMask = (1UL << windowBitCount) - 1UL; return exponent.Low & directMask; } - } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs index 59c7e42..e86e30b 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs @@ -147,13 +147,11 @@ public class MersenneNumberLucasLehmerGpuTester return false; } - GpuPrimeWorkLimiter.Acquire(); bool result; if (exponent >= 128UL) { if (!TryGetNttParameters(exponent, out var nttMod, out var primitiveRoot)) { - GpuPrimeWorkLimiter.Release(); throw new NotSupportedException("NTT parameters for the given exponent are not available."); } @@ -161,6 +159,7 @@ public class MersenneNumberLucasLehmerGpuTester } else { + GpuPrimeWorkLimiter.Acquire(); var accelerator = AcceleratorPool.Shared.Rent(); var stream = accelerator.CreateStream(); var modulus = new GpuUInt128(((UInt128)1 << (int)exponent) - 1UL); // TODO: Cache these Mersenne moduli per exponent so LL GPU runs skip rebuilding them every launch. diff --git a/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs b/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs index e5736b2..7463757 100644 --- a/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs +++ b/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs @@ -10,25 +10,28 @@ public static class Pow2MontgomeryGpuCalculator [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong CalculateKeep(ulong exponent, in MontgomeryDivisorData divisor) { + GpuPrimeWorkLimiter.Acquire(); var gpu = Pow2MontgomeryAccelerator.Rent(1); Accelerator? accelerator = gpu.Accelerator; - AcceleratorStream stream = gpu.Stream!; + AcceleratorStream stream = AcceleratorStreamPool.Rent(accelerator); + // AcceleratorStream stream = gpu.Stream!; MemoryBuffer1D exponentBuffer = gpu.Input; - MemoryBuffer1D resultBuffer = gpu.Output; exponentBuffer.View.CopyFromCPU(stream, ref exponent, 1); // We don't need to worry about any left-overs here. // resultBuffer.MemSetToZero(stream); + MemoryBuffer1D resultBuffer = gpu.OutputUlong; var kernel = gpu.KeepMontgomeryKernel; + var kernelLauncher = kernel.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); - GpuPrimeWorkLimiter.Acquire(); - kernel(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); + kernelLauncher(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); ulong result = 0UL; resultBuffer.View.CopyToCPU(stream, ref result, 1); stream.Synchronize(); + AcceleratorStreamPool.Return(stream); GpuPrimeWorkLimiter.Release(); Pow2MontgomeryAccelerator.Return(gpu); @@ -38,24 +41,26 @@ public static class Pow2MontgomeryGpuCalculator [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong CalculateConvert(ulong exponent, in MontgomeryDivisorData divisor) { + GpuPrimeWorkLimiter.Acquire(); var gpu = Pow2MontgomeryAccelerator.Rent(1); - AcceleratorStream stream = gpu.Stream!; + AcceleratorStream stream = AcceleratorStreamPool.Rent(gpu.Accelerator); + // AcceleratorStream stream = gpu.Stream!; MemoryBuffer1D exponentBuffer = gpu.Input; - MemoryBuffer1D resultBuffer = gpu.Output; exponentBuffer.View.CopyFromCPU(stream, ref exponent, 1); // We don't need to worry about any left-overs here. // resultBuffer.MemSetToZero(stream); + MemoryBuffer1D resultBuffer = gpu.OutputUlong; var kernel = gpu.ConvertToStandardKernel; - GpuPrimeWorkLimiter.Acquire(); - kernel(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); + kernel.Launch(stream, 1, 1, exponentBuffer.View, divisor, resultBuffer.View); ulong result = 0UL; resultBuffer.View.CopyToCPU(stream, ref result, 1); stream.Synchronize(); + AcceleratorStreamPool.Return(stream); GpuPrimeWorkLimiter.Release(); Pow2MontgomeryAccelerator.Return(gpu); @@ -65,20 +70,24 @@ public static class Pow2MontgomeryGpuCalculator [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong CalculateConvert(Pow2MontgomeryAccelerator gpu, ulong exponent, in MontgomeryDivisorData divisor) { - AcceleratorStream stream = gpu.Stream!; + AcceleratorStream stream = AcceleratorStreamPool.Rent(gpu.Accelerator); + // AcceleratorStream stream = gpu.Stream!; MemoryBuffer1D exponentBuffer = gpu.Input; - MemoryBuffer1D resultBuffer = gpu.Output; exponentBuffer.View.CopyFromCPU(stream, ref exponent, 1); // We don't need to worry about any left-overs here. // resultBuffer.MemSetToZero(stream); + MemoryBuffer1D resultBuffer = gpu.OutputUlong; var kernel = gpu.ConvertToStandardKernel; - kernel(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); + var kernelLauncher = kernel.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); + + kernelLauncher(stream, 1,exponentBuffer.View, divisor, resultBuffer.View); ulong result = 0UL; resultBuffer.View.CopyToCPU(stream, ref result, 1); stream.Synchronize(); + AcceleratorStreamPool.Return(stream); return result; } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index d53db13..481f22e 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -224,6 +224,7 @@ public sealed class HeuristicCombinedPrimeTester private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10) { + GpuPrimeWorkLimiter.Acquire(); var gpu = HeuristicPrimeTesterAccelerator.Rent(1); var accelerator = gpu.Accelerator; var stream = gpu.Stream; @@ -238,7 +239,6 @@ public sealed class HeuristicCombinedPrimeTester flagView1D.CopyFromCPU(stream, ref compositeFlag, 1); stream!.Synchronize(); - GpuPrimeWorkLimiter.Acquire(); kernel( stream, diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index 35447fc..19ae607 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -4,6 +4,7 @@ using System.Numerics; using UInt128 = System.UInt128; using ILGPU; using System.Buffers; +using PerfectNumbers.Core.Cpu; namespace PerfectNumbers.Core; @@ -362,8 +363,8 @@ public static class MersennePrimeFactorTester return 0UL; } - c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; + x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; y = x; d = 1UL; @@ -505,8 +506,8 @@ public static class MersennePrimeFactorTester return 0; } - c = (DeterministicRandom.NextUInt128() % (n - 1)) + 1; - x = (DeterministicRandom.NextUInt128() % (n - 2)) + 2; + c = (DeterministicRandomCpu.NextUInt128() % (n - 1)) + 1; + x = (DeterministicRandomCpu.NextUInt128() % (n - 2)) + 2; y = x; d = 1; diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 9c89384..bc99e91 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -12,7 +12,8 @@ public static class PerfectNumberConstants public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxOddPowersCount = 128; public const int PooledArrayThreshold = 64; - public static readonly int RollingAccelerators = 128; //SharedGpuContext.Device.MaxNumThreadsPerGroup; + public static readonly int RollingAccelerators = 256; //SharedGpuContext.Device.MaxNumThreadsPerGroup; + public const int ThreadsByStream = 1; // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal // divisor-cycle datasets we generate offline. // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index e9b141e..5c16d14 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -3,6 +3,7 @@ using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; +using PerfectNumbers.Core.Cpu; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; @@ -1162,7 +1163,8 @@ internal static partial class PrimeOrderCalculator continue; } - if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) + if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) + // if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) { pending.Add(new PendingEntry(composite, knownComposite: true)); continue; @@ -1492,7 +1494,7 @@ internal static partial class PrimeOrderCalculator return baseTicks + (remainder / 1000L); } - private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor) + private static bool TryPollardRhoCpu(ulong n, long deadlineTimestamp, out ulong factor) { factor = 0UL; if ((n & 1UL) == 0UL) @@ -1510,8 +1512,8 @@ internal static partial class PrimeOrderCalculator return false; } - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; ulong y = x; ulong d = 1UL; @@ -1523,9 +1525,9 @@ internal static partial class PrimeOrderCalculator return false; } - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); + x = AdvancePolynomialCpu(x, c, n); + y = AdvancePolynomialCpu(y, c, n); + y = AdvancePolynomialCpu(y, c, n); ulong diff = x > y ? x - y : y - x; d = BinaryGcd(diff, n); } @@ -1539,7 +1541,7 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) + private static ulong AdvancePolynomialCpu(ulong x, ulong c, ulong modulus) { UInt128 value = (UInt128)x * x + c; return (ulong)(value % modulus); @@ -1810,16 +1812,16 @@ internal static partial class PrimeOrderCalculator while (true) { - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; ulong y = x; ulong d = 1UL; while (d == 1UL) { - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); + x = AdvancePolynomialCpu(x, c, n); + y = AdvancePolynomialCpu(y, c, n); + y = AdvancePolynomialCpu(y, c, n); ulong diff = x > y ? x - y : y - x; d = BinaryGcd(diff, n); } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index a066b06..89a49e0 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,3 +1,5 @@ +using System.Diagnostics; +using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; @@ -37,19 +39,21 @@ internal static partial class PrimeOrderCalculator entries.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); - AcceleratorStream stream = gpu.Stream!; + // AcceleratorStream stream = gpu.Stream!; gpu.EnsureCapacity(entryCount, 1); + + GpuPrimeWorkLimiter.Acquire(); + AcceleratorStream stream = AcceleratorStreamPool.Rent(gpu.Accelerator); gpu.Pow2ModEntriesToTestOnDevice.View.CopyFromCPU(stream, entries); ulong order = phi; var kernel = gpu.CheckFactorsKernel; stream.Synchronize(); - GpuPrimeWorkLimiter.Acquire(); - - kernel.Launch(stream, 1, entryCount, phi, gpu.Pow2ModEntriesToTestOnDevice.View, divisorData, gpu.Output.View); - gpu.Output.View.CopyToCPU(stream, ref order, 1); + kernel.Launch(stream, 1, entryCount, phi, gpu.Pow2ModEntriesToTestOnDevice.View, divisorData, gpu.OutputUlong.View); + gpu.OutputUlong.View.CopyToCPU(stream, ref order, 1); stream.Synchronize(); + AcceleratorStreamPool.Return(stream); GpuPrimeWorkLimiter.Release(); return order; @@ -186,5 +190,34 @@ internal static partial class PrimeOrderCalculator GpuPrimeWorkLimiter.Release(); return result != 0; } + + private static bool TryPollardRhoGpu(Pow2MontgomeryAccelerator gpu, ulong n, out ulong factor) + { + // var stream = gpu.Stream!; + Span randomStateSpan = stackalloc ulong[1]; + randomStateSpan[0] = ThreadStaticDeterministicRandomGpu.Exclusive.State; + + AcceleratorStream stream = AcceleratorStreamPool.Rent(gpu.Accelerator); + gpu.Input.View.CopyFromCPU(stream, randomStateSpan); + + Span factoredSpan = stackalloc byte[1]; + Span factorSpan = stackalloc ulong[1]; + var kernelLauncher = gpu.PollardRhoKernel.CreateLauncherDelegate, ArrayView1D, ArrayView1D>>(); + + kernelLauncher(stream, n, 1, gpu.Input.View, gpu.OutputByte.View, gpu.OutputUlong.View); + + gpu.OutputByte.View.CopyToCPU(stream, factoredSpan); + gpu.OutputUlong.View.CopyToCPU(stream, factorSpan); + gpu.Input.View.CopyToCPU(stream, randomStateSpan); + stream.Synchronize(); + + AcceleratorStreamPool.Return(stream); + ThreadStaticDeterministicRandomGpu.Exclusive.SetState(randomStateSpan[0]); + + bool factored = factoredSpan[0] != 0; + factor = factored ? factorSpan[0] : 0UL; + return factored; + } + } diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index c6dc655..4dd3ff5 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -93,23 +93,21 @@ public sealed class PrimeTester [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool IsPrimeGpu(ulong n) { - GpuPrimeWorkLimiter.Acquire(); var gpu = Pow2MontgomeryAccelerator.Rent(1); - var accelerator = gpu.Accelerator; - var stream = gpu.Stream!; - var kernel = gpu.SmallPrimeSieveKernel; + var kernel = gpu.SmallPrimeSieveKernel!; - ulong value = n; // Span flag = stackalloc byte[1]; - byte flag = 0; - var input = gpu.PrimeTestInput; - var output = gpu.PrimeTestOutput; - var inputView = input.View; - var outputView = output.View; + var inputView = gpu.PrimeTestInput.View; - inputView.CopyFromCPU(stream, ref value, 1); - kernel( + GpuPrimeWorkLimiter.Acquire(); + AcceleratorStream stream = AcceleratorStreamPool.Rent(gpu.Accelerator); + // var stream = gpu.Stream!; + inputView.CopyFromCPU(stream, ref n, 1); + + var outputView = gpu.OutputByte.View; + + kernel.Launch( stream, 1, inputView, @@ -122,9 +120,12 @@ public sealed class PrimeTester gpu.DevicePrimesPow2LastThree.View, gpu.DevicePrimesPow2LastNine.View, outputView); + + byte flag = 0; outputView.CopyToCPU(stream, ref flag, 1); stream.Synchronize(); + AcceleratorStreamPool.Return(stream); // gpu.Dispose(); Pow2MontgomeryAccelerator.Return(gpu); GpuPrimeWorkLimiter.Release(); @@ -135,21 +136,21 @@ public sealed class PrimeTester [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool IsPrimeGpu(Pow2MontgomeryAccelerator gpu, ulong n) { - var accelerator = gpu.Accelerator; - var stream = gpu.Stream!; - Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>? kernel = gpu.SmallPrimeSieveKernel; - - ulong value = n; - // Span flag = stackalloc byte[1]; - byte flag = 0; - + // var stream = gpu.Stream!; var input = gpu.PrimeTestInput; - var output = gpu.PrimeTestOutput; var inputView = input.View; + + AcceleratorStream stream = AcceleratorStreamPool.Rent(gpu.Accelerator); + inputView.CopyFromCPU(stream, ref n, 1); + + // Span flag = stackalloc byte[1]; + var output = gpu.OutputByte; var outputView = output.View; - inputView.CopyFromCPU(stream, ref value, 1); - kernel( + Kernel kernel = gpu.SmallPrimeSieveKernel!; + var kernelLauncher = kernel.CreateLauncherDelegate, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>(); + + kernelLauncher( stream, 1, inputView, @@ -162,8 +163,11 @@ public sealed class PrimeTester gpu.DevicePrimesPow2LastThree.View, gpu.DevicePrimesPow2LastNine.View, outputView); + + byte flag = 0; outputView.CopyToCPU(stream, ref flag, 1); stream.Synchronize(); + AcceleratorStreamPool.Return(stream); return flag != 0; } @@ -199,27 +203,27 @@ public sealed class PrimeTester GpuPrimeWorkLimiter.Acquire(); var gpu = Pow2MontgomeryAccelerator.Rent(GpuBatchSize); var accelerator = gpu.Accelerator; - var stream = gpu.Stream!; - var kernel = gpu.SmallPrimeSieveKernel; + // var stream = gpu.Stream!; + var kernel = gpu.SmallPrimeSieveKernel!; int totalLength = values.Length; int batchSize = GpuBatchSize; var input = gpu.PrimeTestInput; - var output = gpu.PrimeTestOutput; + var output = gpu.OutputByte; var inputView = input.View; var outputView = output.View; int pos = 0; + AcceleratorStream stream = AcceleratorStreamPool.Rent(accelerator); while (pos < totalLength) { int remaining = totalLength - pos; int count = remaining > batchSize ? batchSize : remaining; var valueSlice = values.Slice(pos, count); - ref ulong valueRef = ref MemoryMarshal.GetReference(valueSlice); - inputView.CopyFromCPU(stream, ref valueRef, count); + inputView.CopyFromCPU(stream, valueSlice); - kernel( + kernel.Launch( stream, count, inputView, @@ -241,6 +245,7 @@ public sealed class PrimeTester } stream.Synchronize(); + AcceleratorStreamPool.Return(stream); Pow2MontgomeryAccelerator.Return(gpu); // gpu.Dispose(); GpuPrimeWorkLimiter.Release(); @@ -281,7 +286,7 @@ public sealed class PrimeTester // and divisor-cycle staging instead of allocating new device buffers per call. var gpu = Pow2MontgomeryAccelerator.Rent(1); var accelerator = gpu.Accelerator; - var stream = gpu.Stream!; + // var stream = gpu.Stream!; int length = values.Length; ArrayPool pool = ThreadStaticPools.UlongPool; @@ -297,11 +302,14 @@ public sealed class PrimeTester resultBuffer = accelerator.Allocate1D(length); } + AcceleratorStream stream = AcceleratorStreamPool.Rent(accelerator); inputBuffer.View.CopyFromCPU(stream, ref temp[0], length); var kernel = GetSharesFactorKernel(accelerator); kernel(stream, length, inputBuffer.View, resultBuffer.View); resultBuffer.View.CopyToCPU(stream, in results); stream.Synchronize(); + + AcceleratorStreamPool.Return(stream); pool.Return(temp, clearArray: false); resultBuffer.Dispose(); inputBuffer.Dispose();