diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index 3e2ffd1..7858749 100644 --- a/EvenPerfectBitScanner.Tests/ProgramTests.cs +++ b/EvenPerfectBitScanner.Tests/ProgramTests.cs @@ -73,14 +73,14 @@ public class ProgramTests var useDivisorField = typeof(Program).GetField("_useDivisor", BindingFlags.NonPublic | BindingFlags.Static)!; var divisorField = typeof(Program).GetField("_divisor", BindingFlags.NonPublic | BindingFlags.Static)!; var testerField = typeof(Program).GetField("_divisorTester", BindingFlags.NonPublic | BindingFlags.Static)!; - var candidatesField = typeof(Program).GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)!; + var candidatesField = typeof(MersenneNumberDivisorGpuTester).GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var residueField = typeof(Program).GetField("PResidue", BindingFlags.NonPublic | BindingFlags.Static)!; var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu"); useDivisorField.SetValue(null, true); divisorField.SetValue(null, 0UL); - testerField.SetValue(null, null); + testerField.SetValue(null, new MersenneNumberDivisorGpuTester()); primeField.SetValue(null, new ThreadLocal(() => new PrimeTester(), trackAllValues: true)); residueField.SetValue(null, new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true)); candidatesField.SetValue(null, new (ulong, uint)[] { (7UL, 3U), (23UL, 11U) }); @@ -90,6 +90,7 @@ public class ProgramTests { Program.IsEvenPerfectCandidate(11UL, 32).Should().BeFalse(); Program.IsEvenPerfectCandidate(5UL, 32).Should().BeTrue(); + Program.IsEvenPerfectCandidate(127UL, 32).Should().BeTrue(); } finally { diff --git a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs index 622448f..9b2b6a9 100644 --- a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs @@ -12,6 +12,7 @@ public class PrimeTesterIsPrimeTests [InlineData(2UL, true)] [InlineData(3UL, true)] [InlineData(4UL, false)] + [InlineData(5UL, true)] [InlineData(7UL, true)] [InlineData(97UL, true)] [InlineData(121UL, false)] diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index d386c6f..d94418a 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -1,4 +1,5 @@ using System.Collections.Concurrent; +using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; using static PerfectNumbers.Core.Gpu.GpuContextPool; @@ -29,7 +30,7 @@ public sealed class KernelContainer public Action, ArrayView1D>? Incremental; public Action, ArrayView1D>? Pow2Mod; + ResidueAutomatonArgs, ArrayView, ArrayView1D, ArrayView1D, ArrayView1D>? Pow2Mod; public Action, ArrayView1D>? IncrementalOrder; public Action? SmallCycles; + public MemoryBuffer1D? SmallPrimes; + public MemoryBuffer1D? SmallPrimesPow2; public static T InitOnce(ref T? slot, Func factory) where T : class { @@ -77,17 +80,17 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } } - public readonly Action, ArrayView1D> Pow2ModKernel - { - get - { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModKernelScan)); - } - } - } + public readonly Action, ArrayView1D, ArrayView1D, ArrayView1D> Pow2ModKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan)); + } + } + } public readonly Action, ArrayView1D> IncrementalKernel { @@ -274,49 +277,45 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, - ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView orders, - ArrayView1D smallCycles) + ResidueAutomatonArgs ra, ArrayView orders, + ArrayView1D smallCycles, ArrayView1D smallPrimes, + ArrayView1D smallPrimesPow2) { ulong idx = (ulong)index.X; - ulong step10 = ((exponent % 10UL) << 1) % 10UL; - ulong r10 = q0m10 + (step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong step8 = ((exponent & 7UL) << 1) & 7UL; - ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } + ulong r10 = ra.Q0M10 + (ra.Step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } else { - ulong step3 = ((exponent % 3UL) << 1) % 3UL; - ulong r3 = q0m3 + (step3 * (idx % 3UL)); r3 -= (r3 >= 3UL) ? 3UL : 0UL; + ulong r3 = ra.Q0M3 + (ra.Step3 * (idx % 3UL)); r3 -= (r3 >= 3UL) ? 3UL : 0UL; if (r3 == 0UL) { shouldCheck = false; } else { - ulong step5Loc = ((exponent % 5UL) << 1) % 5UL; - ulong r5 = q0m5 + (step5Loc * (idx % 5UL)); if (r5 >= 5UL) r5 -= 5UL; + ulong r5 = ra.Q0M5 + (ra.Step5 * (idx % 5UL)); if (r5 >= 5UL) r5 -= 5UL; if (r5 == 0UL) { shouldCheck = false; } } } - } - if (!shouldCheck) - { - orders[index] = 0UL; - return; - } + } + if (!shouldCheck) + { + orders[index] = 0UL; + return; + } GpuUInt128 k = kStart + (GpuUInt128)idx; GpuUInt128 q = twoP.Mul64(k) + GpuUInt128.One; - // Small-cycles in-kernel early rejection from device table if (q.High == 0UL && q.Low <= (ulong)(smallCycles.Length - 1)) { uint cyc = smallCycles[(int)q.Low]; @@ -332,9 +331,54 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } + int primesLen = (int)smallPrimes.Length; + for (int i = 0; i < primesLen; i++) + { + ulong square = smallPrimesPow2[i]; + if (new GpuUInt128(0UL, square) > q) + { + break; + } + ulong prime = smallPrimes[i]; + if (Mod128By64(q, prime) == 0UL) + { + orders[index] = 0UL; + return; + } + } + orders[index] = exponent; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Mod128By64(GpuUInt128 value, ulong modulus) + { + ulong rem = 0UL; + ulong part = value.High; + for (int i = 0; i < 64; i++) + { + rem = (rem << 1) | (part >> 63); + if (rem >= modulus) + { + rem -= modulus; + } + part <<= 1; + } + + part = value.Low; + for (int i = 0; i < 64; i++) + { + rem = (rem << 1) | (part >> 63); + if (rem >= modulus) + { + rem -= modulus; + } + part <<= 1; + } + + return rem; + } + // Device-friendly small cycle length for 2 mod divisor private static ulong CalculateCycleLengthSmall(ulong divisor) { @@ -527,6 +571,33 @@ public class GpuKernelPool } } + public static (ArrayView1D primes, ArrayView1D primesPow2) EnsureSmallPrimesOnDevice(Accelerator accelerator) + { + var kernels = GetKernels(accelerator); + if (kernels.SmallPrimes is { } p && kernels.SmallPrimesPow2 is { } p2) + { + return (p.View, p2.View); + } + + lock (kernels) + { + if (kernels.SmallPrimes is { } existing && kernels.SmallPrimesPow2 is { } existingPow2) + { + return (existing.View, existingPow2.View); + } + + var hostPrimes = PrimesGenerator.SmallPrimes; + var hostPow2 = PrimesGenerator.SmallPrimesPow2; + var devicePrimes = accelerator.Allocate1D(hostPrimes.Length); + devicePrimes.View.CopyFromCPU(hostPrimes); + var devicePow2 = accelerator.Allocate1D(hostPow2.Length); + devicePow2.View.CopyFromCPU(hostPow2); + kernels.SmallPrimes = devicePrimes; + kernels.SmallPrimesPow2 = devicePow2; + return (devicePrimes.View, devicePow2.View); + } + } + public static GpuKernelLease GetKernel(bool useGpuOrder) { var limiter = GpuPrimeWorkLimiter.Acquire(); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index 7a8a52a..13819e8 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -7,10 +7,11 @@ namespace PerfectNumbers.Core.Gpu; public sealed class MersenneNumberDivisorGpuTester { - private readonly ConcurrentDictionary>> _kernelCache = new(); + private readonly ConcurrentDictionary>> _kernelCache = new(); + private readonly ConcurrentDictionary> _resultBuffers = new(); - private Action> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(Kernel)); + private Action> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(Kernel)); public static void BuildDivisorCandidates() { @@ -36,15 +37,14 @@ public sealed class MersenneNumberDivisorGpuTester { var gpu = GpuContextPool.RentPreferred(preferCpu: false); var accelerator = gpu.Accelerator; - var kernel = GetKernel(accelerator); - var resultBuffer = accelerator.Allocate1D(1); - kernel(1, exponent, divisor, resultBuffer.View); - accelerator.Synchronize(); - bool divisible = resultBuffer.GetAsArray1D()[0] != 0; - resultBuffer.Dispose(); - gpu.Dispose(); - return divisible; - } + var kernel = GetKernel(accelerator); + var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D(1)); + kernel(1, exponent, divisor, resultBuffer.View); + accelerator.Synchronize(); + bool divisible = resultBuffer.GetAsArray1D()[0] != 0; + gpu.Dispose(); + return divisible; + } private static (ulong divisor, uint cycle)[]? _divisorCandidates; @@ -92,9 +92,13 @@ public sealed class MersenneNumberDivisorGpuTester break; } - d = kMul2 * p + 1UL; - // kkMul2 now becomes cycle. We're reusing existing variable for the best performance - kMul2 = divisorCycles.GetCycle(d); + d = kMul2 * p + 1UL; + if (p < 64UL && d == ((1UL << (int)p) - 1UL)) + { + continue; + } + // kkMul2 now becomes cycle. We're reusing existing variable for the best performance + kMul2 = divisorCycles.GetCycle(d); if (p % kMul2 != 0UL) { continue; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs index 743b0e5..86c2e9e 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs @@ -1,5 +1,6 @@ using System.Buffers; using System.Runtime.InteropServices; +using ILGPU; using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu; @@ -18,9 +19,20 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL; byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters - var kernel = (_kernelType == GpuKernelType.Pow2Mod) - ? gpuLease.Pow2ModKernel - : gpuLease.IncrementalKernel; + var pow2Kernel = gpuLease.Pow2ModKernel; + var incKernel = gpuLease.IncrementalKernel; + ulong step10 = ((exponent % 10UL) << 1) % 10UL; + ulong step8 = ((exponent & 7UL) << 1) & 7UL; + ulong step3 = ((exponent % 3UL) << 1) % 3UL; + ulong step5 = ((exponent % 5UL) << 1) % 5UL; + GpuUInt128 twoPGpu = (GpuUInt128)twoP; + var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); + ArrayView1D smallPrimesView = default; + ArrayView1D smallPrimesPow2View = default; + if (_kernelType == GpuKernelType.Pow2Mod) + { + (smallPrimesView, smallPrimesPow2View) = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); + } var orderBuffer = accelerator.Allocate1D(batchSize); // Avoid giant stack allocations that can trigger StackOverflow when batchSize is large. @@ -39,27 +51,39 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u Span orders = orderArray.AsSpan(0, currentSize); // Precompute residue automaton bases for this batch UInt128 q0 = twoP * kStart + UInt128.One; - ulong q0m10 = q0.Mod10(); - ulong q0m8 = q0.Mod8(); - ulong q0m3 = q0.Mod3(); - ulong q0m5 = q0.Mod5(); - ulong step10 = ((exponent % 10UL) << 1) % 10UL; - ulong step8 = ((exponent & 7UL) << 1) & 7UL; - ulong step3 = ((exponent % 3UL) << 1) % 3UL; - ulong step5 = ((exponent % 5UL) << 1) % 5UL; - - var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); - kernel(currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul, - q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView); + if (_kernelType == GpuKernelType.Pow2Mod) + { + q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); + var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5); + pow2Kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul, + ra, orderBuffer.View, smallCyclesView, smallPrimesView, smallPrimesPow2View); + } + else + { + q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); + incKernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul, + q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView); + } accelerator.Synchronize(); orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); - q = twoP * kStart + UInt128.One; - for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP) + if (_kernelType == GpuKernelType.Pow2Mod) { - if (orders[i] != 0UL) + for (i = 0; i < currentSize; i++) { - if (q.IsPrimeCandidate()) + if (orders[i] != 0UL) + { + Volatile.Write(ref isPrime, false); + break; + } + } + } + else + { + q = twoP * kStart + UInt128.One; + for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP) + { + if (orders[i] != 0UL && q.IsPrimeCandidate()) { Volatile.Write(ref isPrime, false); } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index f836264..5a69cb4 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -15,26 +15,32 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) { var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); - var accelerator = gpuLease.Accelerator; - // Ensure device has small cycles table for in-kernel lookup - var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); - int batchSize = GpuConstants.ScanBatchSize; - UInt128 kStart = 1UL; - byte last = lastIsSeven ? (byte)1 : (byte)0; - var kernel = gpuLease.Pow2ModKernel; + var accelerator = gpuLease.Accelerator; + // Ensure device has small cycles and primes tables for in-kernel lookup + var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); + var (smallPrimesView, smallPrimesPow2View) = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); + int batchSize = GpuConstants.ScanBatchSize; + UInt128 kStart = 1UL; + byte last = lastIsSeven ? (byte)1 : (byte)0; + var kernel = gpuLease.Pow2ModKernel; + ulong step10 = ((exponent % 10UL) << 1) % 10UL; + ulong step8 = ((exponent & 7UL) << 1) & 7UL; + ulong step3 = ((exponent % 3UL) << 1) % 3UL; + ulong step5 = ((exponent % 5UL) << 1) % 5UL; + GpuUInt128 twoPGpu = (GpuUInt128)twoP; - var orderBuffer = accelerator.Allocate1D(batchSize); - ulong[] orderArray = ArrayPool.Shared.Rent(batchSize); - UInt128 remaining; - int currentSize; - int i; - UInt128 batchSize128 = (UInt128)batchSize, q; - Span orders = orderArray.AsSpan(0, batchSize); - try - { - while (kStart < maxK && Volatile.Read(ref isPrime)) - { - remaining = maxK - kStart; + var orderBuffer = accelerator.Allocate1D(batchSize); + ulong[] orderArray = ArrayPool.Shared.Rent(batchSize); + UInt128 remaining; + int currentSize; + int i; + UInt128 batchSize128 = (UInt128)batchSize, q; + Span orders = orderArray.AsSpan(0, batchSize); + try + { + while (kStart < maxK && Volatile.Read(ref isPrime)) + { + remaining = maxK - kStart; if (remaining > batchSize128) { currentSize = batchSize; @@ -45,39 +51,35 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) orders = orderArray.AsSpan(0, currentSize); } - q = twoP * kStart + UInt128.One; - q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); - ulong step10 = ((exponent % 10UL) << 1) % 10UL; - ulong step8 = ((exponent & 7UL) << 1) & 7UL; - ulong step3 = ((exponent % 3UL) << 1) % 3UL; - ulong step5 = ((exponent % 5UL) << 1) % 5UL; + q = twoP * kStart + UInt128.One; + q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); + var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5); - kernel(currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, 0UL, - q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView); + kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL, + ra, orderBuffer.View, smallCyclesView, smallPrimesView, smallPrimesPow2View); - accelerator.Synchronize(); - orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); - accelerator.Synchronize(); - if (!Volatile.Read(ref isPrime)) - { - break; - } + accelerator.Synchronize(); + orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); + if (!Volatile.Read(ref isPrime)) + { + break; + } - for (i = 0; i < currentSize; i++, q += twoP) - { - if (orders[i] != 0UL && q.IsPrimeCandidate()) - { - Volatile.Write(ref isPrime, false); - break; - } - } + for (i = 0; i < currentSize; i++) + { + if (orders[i] != 0UL) + { + Volatile.Write(ref isPrime, false); + break; + } + } - kStart += batchSize128; - } - } - finally - { - ArrayPool.Shared.Return(orderArray); + kStart += batchSize128; + } + } + finally + { + ArrayPool.Shared.Return(orderArray); orderBuffer.Dispose(); gpuLease.Dispose(); } diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index eebeddb..d627bb3 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -42,7 +42,7 @@ public sealed class PrimeTester(bool useInternal = false) { result = false; } - else if (n % 5 == 0) + else if (n > 5UL && (n % 5UL) == 0UL) { result = false; }