diff --git a/PerfectNumbers.Core.Tests/Gpu/LucasLehmerPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/Gpu/LucasLehmerPrimeTesterTests.cs index 79ac24c..7fd6cf1 100644 --- a/PerfectNumbers.Core.Tests/Gpu/LucasLehmerPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/LucasLehmerPrimeTesterTests.cs @@ -29,6 +29,18 @@ public class LucasLehmerPrimeTesterTests tester.IsMersennePrime(exponent).Should().BeFalse(); } + [Fact] + [Trait("Category", "Fast")] + public void IsMersennePrimeBatch_matches_individual() + { + ulong[] exponents = new ulong[] { 3UL, 5UL, 7UL, 11UL, 13UL, 17UL, 19UL, 23UL }; + bool[] expected = new bool[] { true, true, true, false, true, true, true, false }; + var tester = new MersenneNumberLucasLehmerGpuTester(); + bool[] results = new bool[exponents.Length]; + tester.IsMersennePrimeBatch(exponents, results); + results.Should().Equal(expected); + } + [Theory] [InlineData(132UL)] [InlineData(1000UL)] diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs index 96719bc..efcbc8b 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs @@ -1,7 +1,9 @@ +using System.Buffers; using System.Collections.Concurrent; +using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; -using System.Runtime.CompilerServices; +using UInt128 = System.UInt128; namespace PerfectNumbers.Core.Gpu; @@ -12,6 +14,7 @@ public class MersenneNumberLucasLehmerGpuTester private readonly ConcurrentDictionary>> SubSmallKernelCache = new(); private readonly ConcurrentDictionary>> ReduceKernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> IsZeroKernelCache = new(); + private readonly ConcurrentDictionary, ArrayView, ArrayView>> BatchKernelCache = new(); private Action> GetKernel(Accelerator accelerator) { @@ -30,6 +33,9 @@ public class MersenneNumberLucasLehmerGpuTester private Action, ArrayView> GetIsZeroKernel(Accelerator accelerator) => IsZeroKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(IsZeroKernel)); + private Action, ArrayView, ArrayView> GetBatchKernel(Accelerator accelerator) => + BatchKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(KernelBatch)); + // Configurable LL slice size to keep kernels short. Default 32. public int SliceSize = 32; @@ -85,6 +91,76 @@ public class MersenneNumberLucasLehmerGpuTester return result; } + public void ComputeResidues(ReadOnlySpan exponents, Span residues) + { + int count = exponents.Length; + if (residues.Length < count) + { + throw new ArgumentException("Result span too small", nameof(residues)); + } + + for (int i = 0; i < count; i++) + { + if (exponents[i] >= 128UL) + { + throw new NotSupportedException("Batch residue calculation supports exponents < 128 only."); + } + } + + var gpu = GpuContextPool.RentPreferred(preferCpu: false); + var accelerator = gpu.Accelerator; + var kernel = GetBatchKernel(accelerator); + + var expBuffer = accelerator.Allocate1D(count); + ulong[] expArray = ArrayPool.Shared.Rent(count); + exponents.CopyTo(expArray); + expBuffer.View.CopyFromCPU(ref expArray[0], count); + + GpuUInt128[] modulusArray = ArrayPool.Shared.Rent(count); + for (int i = 0; i < count; i++) + { + modulusArray[i] = new GpuUInt128(((UInt128)1 << (int)exponents[i]) - 1UL); + } + + var modBuffer = accelerator.Allocate1D(count); + modBuffer.View.CopyFromCPU(ref modulusArray[0], count); + + var stateBuffer = accelerator.Allocate1D(count); + kernel(count, expBuffer.View, modBuffer.View, stateBuffer.View); + accelerator.Synchronize(); + stateBuffer.View.CopyToCPU(ref residues[0], count); + + stateBuffer.Dispose(); + modBuffer.Dispose(); + expBuffer.Dispose(); + ArrayPool.Shared.Return(modulusArray); + ArrayPool.Shared.Return(expArray); + gpu.Dispose(); + } + + public void IsMersennePrimeBatch(ReadOnlySpan exponents, Span results) + { + int count = exponents.Length; + if (results.Length < count) + { + throw new ArgumentException("Result span too small", nameof(results)); + } + + GpuUInt128[] buffer = ArrayPool.Shared.Rent(count); + try + { + ComputeResidues(exponents, buffer.AsSpan(0, count)); + for (int i = 0; i < count; i++) + { + results[i] = buffer[i].IsZero; + } + } + finally + { + ArrayPool.Shared.Return(buffer); + } + } + private bool IsMersennePrimeNtt(ulong exponent, GpuUInt128 nttMod, GpuUInt128 primitiveRoot, bool runOnGpu) { var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu); @@ -731,6 +807,21 @@ public class MersenneNumberLucasLehmerGpuTester } } + private static void KernelBatch(Index1D index, ArrayView exponents, ArrayView moduli, ArrayView states) + { + int idx = index.X; + ulong exponent = exponents[idx]; + GpuUInt128 modulus = moduli[idx]; + var s = new GpuUInt128(0UL, 4UL); + ulong limit = exponent - 2UL; + for (ulong i = 0UL; i < limit; i++) + { + s = SquareModMersenne128(s, exponent).SubMod(2UL, modulus); + } + + states[idx] = s; + } + private static void Kernel(Index1D index, ulong exponent, GpuUInt128 modulus, ArrayView state) { // Lucas–Lehmer iteration in the field GF(2^p-1).