diff --git a/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs b/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs index 02fc664..d39c3d1 100644 --- a/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs +++ b/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs @@ -1,4 +1,5 @@ using FluentAssertions; +using PerfectNumbers.Core; using Xunit; namespace PerfectNumbers.Core.Tests; @@ -19,6 +20,45 @@ public class PrimesGeneratorTests { squares[i].Should().Be(primes[i] * primes[i]); } + + primes.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); + squares.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); + } + + [Fact] + [Trait("Category", "Fast")] + public void Residue_prime_tables_follow_digit_rules() + { + var lastOne = PrimesGenerator.SmallPrimesLastOne; + var lastOnePow2 = PrimesGenerator.SmallPrimesPow2LastOne; + var lastSeven = PrimesGenerator.SmallPrimesLastSeven; + var lastSevenPow2 = PrimesGenerator.SmallPrimesPow2LastSeven; + + lastOne.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); + lastSeven.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); + lastOnePow2.Length.Should().Be(lastOne.Length); + lastSevenPow2.Length.Should().Be(lastSeven.Length); + + lastOne[..5].Should().Equal([3U, 7U, 11U, 13U, 19U]); + lastSeven[..5].Should().Equal([3U, 7U, 11U, 13U, 17U]); + + for (int i = 0; i < lastOne.Length; i++) + { + uint prime = lastOne[i]; + uint mod10 = prime % 10U; + (mod10 == 1U || mod10 == 3U || mod10 == 9U || prime == 7U || prime == 11U) + .Should().BeTrue(); + lastOnePow2[i].Should().Be(prime * (ulong)prime); + } + + for (int i = 0; i < lastSeven.Length; i++) + { + uint prime = lastSeven[i]; + uint mod10 = prime % 10U; + (mod10 == 3U || mod10 == 7U || mod10 == 9U || prime == 11U) + .Should().BeTrue(); + lastSevenPow2[i].Should().Be(prime * (ulong)prime); + } } } diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index d94418a..8ef52b0 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -30,7 +30,9 @@ public sealed class KernelContainer public Action, ArrayView1D>? Incremental; public Action, ArrayView1D, ArrayView1D, ArrayView1D>? Pow2Mod; + ResidueAutomatonArgs, ArrayView, ArrayView1D, + ArrayView1D, ArrayView1D, + ArrayView1D, ArrayView1D>? Pow2Mod; public Action, ArrayView1D>? IncrementalOrder; public Action? SmallCycles; - public MemoryBuffer1D? SmallPrimes; - public MemoryBuffer1D? SmallPrimesPow2; + public MemoryBuffer1D? SmallPrimesLastOne; + public MemoryBuffer1D? SmallPrimesPow2LastOne; + public MemoryBuffer1D? SmallPrimesLastSeven; + public MemoryBuffer1D? SmallPrimesPow2LastSeven; - public static T InitOnce(ref T? slot, Func factory) where T : class - { + public static T InitOnce(ref T? slot, Func factory) where T : class + { var current = Volatile.Read(ref slot); if (current is not null) { @@ -55,8 +59,20 @@ public sealed class KernelContainer current = factory(); Volatile.Write(ref slot, current); } - return current; - } + return current; + } +} + +public readonly struct ResiduePrimeViews( + ArrayView1D lastOne, + ArrayView1D lastSeven, + ArrayView1D lastOnePow2, + ArrayView1D lastSevenPow2) +{ + public readonly ArrayView1D LastOne = lastOne; + public readonly ArrayView1D LastSeven = lastSeven; + public readonly ArrayView1D LastOnePow2 = lastOnePow2; + public readonly ArrayView1D LastSevenPow2 = lastSevenPow2; } public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels) @@ -80,14 +96,14 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } } - public readonly Action, ArrayView1D, ArrayView1D, ArrayView1D> Pow2ModKernel + public readonly Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> Pow2ModKernel { get { var accel = Accelerator; lock (_gpu!.Value.KernelInitLock) { - return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan)); + return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan)); } } } @@ -278,8 +294,11 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, ResidueAutomatonArgs ra, ArrayView orders, - ArrayView1D smallCycles, ArrayView1D smallPrimes, - ArrayView1D smallPrimesPow2) + ArrayView1D smallCycles, + ArrayView1D smallPrimesLastOne, + ArrayView1D smallPrimesLastSeven, + ArrayView1D smallPrimesPow2LastOne, + ArrayView1D smallPrimesPow2LastSeven) { ulong idx = (ulong)index.X; ulong r10 = ra.Q0M10 + (ra.Step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL; @@ -331,15 +350,17 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } - int primesLen = (int)smallPrimes.Length; + ArrayView1D primes = lastIsSeven != 0 ? smallPrimesLastSeven : smallPrimesLastOne; + ArrayView1D primesPow2 = lastIsSeven != 0 ? smallPrimesPow2LastSeven : smallPrimesPow2LastOne; + int primesLen = (int)primes.Length; for (int i = 0; i < primesLen; i++) { - ulong square = smallPrimesPow2[i]; + ulong square = primesPow2[i]; if (new GpuUInt128(0UL, square) > q) { break; } - ulong prime = smallPrimes[i]; + ulong prime = primes[i]; if (Mod128By64(q, prime) == 0UL) { orders[index] = 0UL; @@ -571,30 +592,47 @@ public class GpuKernelPool } } - public static (ArrayView1D primes, ArrayView1D primesPow2) EnsureSmallPrimesOnDevice(Accelerator accelerator) + public static ResiduePrimeViews EnsureSmallPrimesOnDevice(Accelerator accelerator) { var kernels = GetKernels(accelerator); - if (kernels.SmallPrimes is { } p && kernels.SmallPrimesPow2 is { } p2) + if (kernels.SmallPrimesLastOne is { } lastOne && + kernels.SmallPrimesLastSeven is { } lastSeven && + kernels.SmallPrimesPow2LastOne is { } lastOnePow2 && + kernels.SmallPrimesPow2LastSeven is { } lastSevenPow2) { - return (p.View, p2.View); + return new ResiduePrimeViews(lastOne.View, lastSeven.View, lastOnePow2.View, lastSevenPow2.View); } lock (kernels) { - if (kernels.SmallPrimes is { } existing && kernels.SmallPrimesPow2 is { } existingPow2) + if (kernels.SmallPrimesLastOne is { } existingLastOne && + kernels.SmallPrimesLastSeven is { } existingLastSeven && + kernels.SmallPrimesPow2LastOne is { } existingLastOnePow2 && + kernels.SmallPrimesPow2LastSeven is { } existingLastSevenPow2) { - return (existing.View, existingPow2.View); + return new ResiduePrimeViews(existingLastOne.View, existingLastSeven.View, existingLastOnePow2.View, existingLastSevenPow2.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); + var hostLastOne = PrimesGenerator.SmallPrimesLastOne; + var hostLastSeven = PrimesGenerator.SmallPrimesLastSeven; + var hostLastOnePow2 = PrimesGenerator.SmallPrimesPow2LastOne; + var hostLastSevenPow2 = PrimesGenerator.SmallPrimesPow2LastSeven; + + var deviceLastOne = accelerator.Allocate1D(hostLastOne.Length); + deviceLastOne.View.CopyFromCPU(hostLastOne); + var deviceLastSeven = accelerator.Allocate1D(hostLastSeven.Length); + deviceLastSeven.View.CopyFromCPU(hostLastSeven); + var deviceLastOnePow2 = accelerator.Allocate1D(hostLastOnePow2.Length); + deviceLastOnePow2.View.CopyFromCPU(hostLastOnePow2); + var deviceLastSevenPow2 = accelerator.Allocate1D(hostLastSevenPow2.Length); + deviceLastSevenPow2.View.CopyFromCPU(hostLastSevenPow2); + + kernels.SmallPrimesLastOne = deviceLastOne; + kernels.SmallPrimesLastSeven = deviceLastSeven; + kernels.SmallPrimesPow2LastOne = deviceLastOnePow2; + kernels.SmallPrimesPow2LastSeven = deviceLastSevenPow2; + + return new ResiduePrimeViews(deviceLastOne.View, deviceLastSeven.View, deviceLastOnePow2.View, deviceLastSevenPow2.View); } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs index 86c2e9e..8dd0352 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs @@ -27,11 +27,10 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u ulong step5 = ((exponent % 5UL) << 1) % 5UL; GpuUInt128 twoPGpu = (GpuUInt128)twoP; var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); - ArrayView1D smallPrimesView = default; - ArrayView1D smallPrimesPow2View = default; + ResiduePrimeViews primeViews = default; if (_kernelType == GpuKernelType.Pow2Mod) { - (smallPrimesView, smallPrimesPow2View) = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); + primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); } var orderBuffer = accelerator.Allocate1D(batchSize); @@ -56,7 +55,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u 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); + ra, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2); } else { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index 5a69cb4..b8cfbd2 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -18,7 +18,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) 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); + ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); int batchSize = GpuConstants.ScanBatchSize; UInt128 kStart = 1UL; byte last = lastIsSeven ? (byte)1 : (byte)0; @@ -56,7 +56,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5); kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL, - ra, orderBuffer.View, smallCyclesView, smallPrimesView, smallPrimesPow2View); + ra, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2); accelerator.Synchronize(); orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); diff --git a/PerfectNumbers.Core/PrimesGenerator.cs b/PerfectNumbers.Core/PrimesGenerator.cs index bf2ab7b..0ec6ed1 100644 --- a/PerfectNumbers.Core/PrimesGenerator.cs +++ b/PerfectNumbers.Core/PrimesGenerator.cs @@ -1,80 +1,113 @@ +using System.Collections.Generic; +using System.Runtime.CompilerServices; + namespace PerfectNumbers.Core; public static class PrimesGenerator { - private static uint[]? _primes; - private static ulong[]? _primesPow2; +public static readonly uint[] SmallPrimes; +public static readonly ulong[] SmallPrimesPow2; +public static readonly uint[] SmallPrimesLastOne; +public static readonly ulong[] SmallPrimesPow2LastOne; +public static readonly uint[] SmallPrimesLastSeven; +public static readonly ulong[] SmallPrimesPow2LastSeven; - public static readonly uint[] SmallPrimes = BuildSmallPrimesInternal(); - public static readonly ulong[] SmallPrimesPow2 = BuildSmallPrimesPow2Internal(); +static PrimesGenerator() +{ +BuildSmallPrimes( +out SmallPrimes, +out SmallPrimesPow2, +out SmallPrimesLastOne, +out SmallPrimesPow2LastOne, +out SmallPrimesLastSeven, +out SmallPrimesPow2LastSeven); +} +private static void BuildSmallPrimes( +out uint[] all, +out ulong[] allPow2, +out uint[] lastOne, +out ulong[] lastOnePow2, +out uint[] lastSeven, +out ulong[] lastSevenPow2) +{ +uint target = PerfectNumberConstants.PrimesLimit; +int targetInt = checked((int)target); +all = new uint[targetInt]; +allPow2 = new ulong[targetInt]; +lastOne = new uint[targetInt]; +lastOnePow2 = new ulong[targetInt]; +lastSeven = new uint[targetInt]; +lastSevenPow2 = new ulong[targetInt]; - private static uint[] BuildSmallPrimesInternal() - { - if (_primes == null) - { - BuildSmallPrimes(); - } +List primes = new(targetInt * 4 / 3 + 16); +uint candidate = 2U; +int allCount = 0; +int lastOneCount = 0; +int lastSevenCount = 0; - return _primes!; - } +while (allCount < targetInt || lastOneCount < targetInt || lastSevenCount < targetInt) +{ +bool isPrime = true; +int primesCount = primes.Count; +for (int i = 0; i < primesCount; i++) +{ +uint p = primes[i]; +if (p * p > candidate) +{ +break; +} - private static ulong[] BuildSmallPrimesPow2Internal() - { - if (_primesPow2 == null) - { - BuildSmallPrimes(); - } +if (candidate % p == 0U) +{ +isPrime = false; +break; +} +} - return _primesPow2!; - } +if (isPrime) +{ +primes.Add(candidate); +if (allCount < targetInt) +{ +all[allCount] = candidate; +allPow2[allCount] = checked((ulong)candidate * candidate); +allCount++; +} - private static void BuildSmallPrimes() - { - uint limit = PerfectNumberConstants.PrimesLimit; - var primes = new uint[limit]; - primes[0] = 2; - var primesPow2 = new ulong[limit]; - primesPow2[0] = 4; +if (candidate != 2U) +{ +if (lastOneCount < targetInt && IsAllowedForLastOne(candidate)) +{ +lastOne[lastOneCount] = candidate; +lastOnePow2[lastOneCount] = checked((ulong)candidate * candidate); +lastOneCount++; +} - int j; - uint i = 3, p, sqrt; - bool isPrime; - int primeIndex = 1; +if (lastSevenCount < targetInt && IsAllowedForLastSeven(candidate)) +{ +lastSeven[lastSevenCount] = candidate; +lastSevenPow2[lastSevenCount] = checked((ulong)candidate * candidate); +lastSevenCount++; +} +} +} - for (; i <= limit; i += 2) - { - isPrime = true; - sqrt = (uint)Math.Sqrt(i); - for (j = 0; j < primeIndex; j++) - { - p = primes[j]; - if (p > sqrt) - { - break; - } +candidate = candidate == 2U ? 3U : candidate + 2U; +} +} - if (i % p == 0UL) - { - isPrime = false; - break; - } - } +[MethodImpl(MethodImplOptions.AggressiveInlining)] +private static bool IsAllowedForLastOne(uint prime) +{ +uint mod10 = prime % 10U; +return mod10 == 1U || mod10 == 3U || mod10 == 9U || prime == 7U || prime == 11U; +} - if (isPrime) - { - primes[primeIndex] = i; - primesPow2[primeIndex] = checked((ulong)i * i); - primeIndex++; - } - } - - Array.Resize(ref primes, primeIndex); - Array.Resize(ref primesPow2, primeIndex); - - _primes = primes; - _primesPow2 = primesPow2; - - return; - } -} \ No newline at end of file +[MethodImpl(MethodImplOptions.AggressiveInlining)] +private static bool IsAllowedForLastSeven(uint prime) +{ +uint mod10 = prime % 10U; +return mod10 == 3U || mod10 == 7U || mod10 == 9U || prime == 11U; +} +}