Refine GPU residue prime tables

This commit is contained in:
Marek Lemańczyk committed 2025-09-17 17:35:26 +02:00
1 parent 88d0c9c271
commit 34cc0b0fdd
5 files changed
+188 -78

No files matched your search

@@ -1,4 +1,5 @@
using FluentAssertions; using FluentAssertions;
using PerfectNumbers.Core;
using Xunit; using Xunit;
namespace PerfectNumbers.Core.Tests; namespace PerfectNumbers.Core.Tests;
@@ -19,6 +20,45 @@ public class PrimesGeneratorTests
{ {
squares[i].Should().Be(primes[i] * primes[i]); 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);
}
} }
} }
+62 -24
View File
@@ -30,7 +30,9 @@ public sealed class KernelContainer
public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>? Incremental; ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>? Incremental;
public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2Mod; ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>,
ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>,
ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2Mod;
public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>? IncrementalOrder; ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>? IncrementalOrder;
public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
@@ -38,8 +40,10 @@ public sealed class KernelContainer
// Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length. // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallCycles; public MemoryBuffer1D<uint, Stride1D.Dense>? SmallCycles;
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimes; public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastOne;
public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2; public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastOne;
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastSeven;
public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastSeven;
public static T InitOnce<T>(ref T? slot, Func<T> factory) where T : class public static T InitOnce<T>(ref T? slot, Func<T> factory) where T : class
{ {
@@ -59,6 +63,18 @@ public sealed class KernelContainer
} }
} }
public readonly struct ResiduePrimeViews(
ArrayView1D<uint, Stride1D.Dense> lastOne,
ArrayView1D<uint, Stride1D.Dense> lastSeven,
ArrayView1D<ulong, Stride1D.Dense> lastOnePow2,
ArrayView1D<ulong, Stride1D.Dense> lastSevenPow2)
{
public readonly ArrayView1D<uint, Stride1D.Dense> LastOne = lastOne;
public readonly ArrayView1D<uint, Stride1D.Dense> LastSeven = lastSeven;
public readonly ArrayView1D<ulong, Stride1D.Dense> LastOnePow2 = lastOnePow2;
public readonly ArrayView1D<ulong, Stride1D.Dense> LastSevenPow2 = lastSevenPow2;
}
public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels) public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels)
{ {
private bool disposedValue; private bool disposedValue;
@@ -80,14 +96,14 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
} }
} }
public readonly Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> Pow2ModKernel public readonly Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> Pow2ModKernel
{ {
get get
{ {
var accel = Accelerator; var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock) lock (_gpu!.Value.KernelInitLock)
{ {
return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernelScan)); return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(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 _, private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _,
ResidueAutomatonArgs ra, ArrayView<ulong> orders, ResidueAutomatonArgs ra, ArrayView<ulong> orders,
ArrayView1D<uint, Stride1D.Dense> smallCycles, ArrayView1D<uint, Stride1D.Dense> smallPrimes, ArrayView1D<uint, Stride1D.Dense> smallCycles,
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2) ArrayView1D<uint, Stride1D.Dense> smallPrimesLastOne,
ArrayView1D<uint, Stride1D.Dense> smallPrimesLastSeven,
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2LastOne,
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2LastSeven)
{ {
ulong idx = (ulong)index.X; ulong idx = (ulong)index.X;
ulong r10 = ra.Q0M10 + (ra.Step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL; 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; return;
} }
int primesLen = (int)smallPrimes.Length; ArrayView1D<uint, Stride1D.Dense> primes = lastIsSeven != 0 ? smallPrimesLastSeven : smallPrimesLastOne;
ArrayView1D<ulong, Stride1D.Dense> primesPow2 = lastIsSeven != 0 ? smallPrimesPow2LastSeven : smallPrimesPow2LastOne;
int primesLen = (int)primes.Length;
for (int i = 0; i < primesLen; i++) for (int i = 0; i < primesLen; i++)
{ {
ulong square = smallPrimesPow2[i]; ulong square = primesPow2[i];
if (new GpuUInt128(0UL, square) > q) if (new GpuUInt128(0UL, square) > q)
{ {
break; break;
} }
ulong prime = smallPrimes[i]; ulong prime = primes[i];
if (Mod128By64(q, prime) == 0UL) if (Mod128By64(q, prime) == 0UL)
{ {
orders[index] = 0UL; orders[index] = 0UL;
@@ -571,30 +592,47 @@ public class GpuKernelPool
} }
} }
public static (ArrayView1D<uint, Stride1D.Dense> primes, ArrayView1D<ulong, Stride1D.Dense> primesPow2) EnsureSmallPrimesOnDevice(Accelerator accelerator) public static ResiduePrimeViews EnsureSmallPrimesOnDevice(Accelerator accelerator)
{ {
var kernels = GetKernels(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) 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 hostLastOne = PrimesGenerator.SmallPrimesLastOne;
var hostPow2 = PrimesGenerator.SmallPrimesPow2; var hostLastSeven = PrimesGenerator.SmallPrimesLastSeven;
var devicePrimes = accelerator.Allocate1D<uint>(hostPrimes.Length); var hostLastOnePow2 = PrimesGenerator.SmallPrimesPow2LastOne;
devicePrimes.View.CopyFromCPU(hostPrimes); var hostLastSevenPow2 = PrimesGenerator.SmallPrimesPow2LastSeven;
var devicePow2 = accelerator.Allocate1D<ulong>(hostPow2.Length);
devicePow2.View.CopyFromCPU(hostPow2); var deviceLastOne = accelerator.Allocate1D<uint>(hostLastOne.Length);
kernels.SmallPrimes = devicePrimes; deviceLastOne.View.CopyFromCPU(hostLastOne);
kernels.SmallPrimesPow2 = devicePow2; var deviceLastSeven = accelerator.Allocate1D<uint>(hostLastSeven.Length);
return (devicePrimes.View, devicePow2.View); deviceLastSeven.View.CopyFromCPU(hostLastSeven);
var deviceLastOnePow2 = accelerator.Allocate1D<ulong>(hostLastOnePow2.Length);
deviceLastOnePow2.View.CopyFromCPU(hostLastOnePow2);
var deviceLastSevenPow2 = accelerator.Allocate1D<ulong>(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);
} }
} }
@@ -27,11 +27,10 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
ulong step5 = ((exponent % 5UL) << 1) % 5UL; ulong step5 = ((exponent % 5UL) << 1) % 5UL;
GpuUInt128 twoPGpu = (GpuUInt128)twoP; GpuUInt128 twoPGpu = (GpuUInt128)twoP;
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ArrayView1D<uint, Stride1D.Dense> smallPrimesView = default; ResiduePrimeViews primeViews = default;
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2View = default;
if (_kernelType == GpuKernelType.Pow2Mod) if (_kernelType == GpuKernelType.Pow2Mod)
{ {
(smallPrimesView, smallPrimesPow2View) = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
} }
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize); var orderBuffer = accelerator.Allocate1D<ulong>(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); 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); var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
pow2Kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul, 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 else
{ {
@@ -18,7 +18,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
var accelerator = gpuLease.Accelerator; var accelerator = gpuLease.Accelerator;
// Ensure device has small cycles and primes tables for in-kernel lookup // Ensure device has small cycles and primes tables for in-kernel lookup
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
var (smallPrimesView, smallPrimesPow2View) = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
int batchSize = GpuConstants.ScanBatchSize; int batchSize = GpuConstants.ScanBatchSize;
UInt128 kStart = 1UL; UInt128 kStart = 1UL;
byte last = lastIsSeven ? (byte)1 : (byte)0; 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); var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL, 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(); accelerator.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
+81 -48
View File
@@ -1,60 +1,64 @@
using System.Collections.Generic;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core; namespace PerfectNumbers.Core;
public static class PrimesGenerator public static class PrimesGenerator
{ {
private static uint[]? _primes; public static readonly uint[] SmallPrimes;
private static ulong[]? _primesPow2; 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(); static PrimesGenerator()
public static readonly ulong[] SmallPrimesPow2 = BuildSmallPrimesPow2Internal();
private static uint[] BuildSmallPrimesInternal()
{ {
if (_primes == null) BuildSmallPrimes(
{ out SmallPrimes,
BuildSmallPrimes(); out SmallPrimesPow2,
out SmallPrimesLastOne,
out SmallPrimesPow2LastOne,
out SmallPrimesLastSeven,
out SmallPrimesPow2LastSeven);
} }
return _primes!; private static void BuildSmallPrimes(
} out uint[] all,
out ulong[] allPow2,
private static ulong[] BuildSmallPrimesPow2Internal() out uint[] lastOne,
out ulong[] lastOnePow2,
out uint[] lastSeven,
out ulong[] lastSevenPow2)
{ {
if (_primesPow2 == null) uint target = PerfectNumberConstants.PrimesLimit;
{ int targetInt = checked((int)target);
BuildSmallPrimes(); all = new uint[targetInt];
} allPow2 = new ulong[targetInt];
lastOne = new uint[targetInt];
lastOnePow2 = new ulong[targetInt];
lastSeven = new uint[targetInt];
lastSevenPow2 = new ulong[targetInt];
return _primesPow2!; List<uint> primes = new(targetInt * 4 / 3 + 16);
} uint candidate = 2U;
int allCount = 0;
int lastOneCount = 0;
int lastSevenCount = 0;
private static void BuildSmallPrimes() while (allCount < targetInt || lastOneCount < targetInt || lastSevenCount < targetInt)
{ {
uint limit = PerfectNumberConstants.PrimesLimit; bool isPrime = true;
var primes = new uint[limit]; int primesCount = primes.Count;
primes[0] = 2; for (int i = 0; i < primesCount; i++)
var primesPow2 = new ulong[limit];
primesPow2[0] = 4;
int j;
uint i = 3, p, sqrt;
bool isPrime;
int primeIndex = 1;
for (; i <= limit; i += 2)
{ {
isPrime = true; uint p = primes[i];
sqrt = (uint)Math.Sqrt(i); if (p * p > candidate)
for (j = 0; j < primeIndex; j++)
{
p = primes[j];
if (p > sqrt)
{ {
break; break;
} }
if (i % p == 0UL) if (candidate % p == 0U)
{ {
isPrime = false; isPrime = false;
break; break;
@@ -63,18 +67,47 @@ public static class PrimesGenerator
if (isPrime) if (isPrime)
{ {
primes[primeIndex] = i; primes.Add(candidate);
primesPow2[primeIndex] = checked((ulong)i * i); if (allCount < targetInt)
primeIndex++; {
all[allCount] = candidate;
allPow2[allCount] = checked((ulong)candidate * candidate);
allCount++;
}
if (candidate != 2U)
{
if (lastOneCount < targetInt && IsAllowedForLastOne(candidate))
{
lastOne[lastOneCount] = candidate;
lastOnePow2[lastOneCount] = checked((ulong)candidate * candidate);
lastOneCount++;
}
if (lastSevenCount < targetInt && IsAllowedForLastSeven(candidate))
{
lastSeven[lastSevenCount] = candidate;
lastSevenPow2[lastSevenCount] = checked((ulong)candidate * candidate);
lastSevenCount++;
}
} }
} }
Array.Resize(ref primes, primeIndex); candidate = candidate == 2U ? 3U : candidate + 2U;
Array.Resize(ref primesPow2, primeIndex); }
}
_primes = primes; [MethodImpl(MethodImplOptions.AggressiveInlining)]
_primesPow2 = primesPow2; private static bool IsAllowedForLastOne(uint prime)
{
uint mod10 = prime % 10U;
return mod10 == 1U || mod10 == 3U || mod10 == 9U || prime == 7U || prime == 11U;
}
return; [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastSeven(uint prime)
{
uint mod10 = prime % 10U;
return mod10 == 3U || mod10 == 7U || mod10 == 9U || prime == 11U;
} }
} }