mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #23 from mlemanczyk/codex/modify-evenperfectbitscanner-for-mersenne-primes

Browse Source 
Refine GPU residue prime tables
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-17 17:36:09 +02:00
2 parents 88d0c9c271 34cc0b0fdd
commit b1c29ad16e
5 files changed
+209 -99
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs
+40
View File
@@ -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);
}
}
}
PerfectNumbers.Core/Gpu/GpuKernelPool.cs
+66 -28
View File
@@ -30,7 +30,9 @@ public sealed class KernelContainer
public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>? Incremental;
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,
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>? IncrementalOrder;
public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
@@ -38,11 +40,13 @@ public sealed class KernelContainer
// Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallCycles;
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimes;
public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2;
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastOne;
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
{
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<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)
@@ -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
{
var accel = Accelerator;
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 _,
ResidueAutomatonArgs ra, ArrayView<ulong> orders,
ArrayView1D<uint, Stride1D.Dense> smallCycles, ArrayView1D<uint, Stride1D.Dense> smallPrimes,
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2)
ArrayView1D<uint, Stride1D.Dense> smallCycles,
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 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<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++)
{
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<uint, Stride1D.Dense> primes, ArrayView1D<ulong, Stride1D.Dense> 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<uint>(hostPrimes.Length);
devicePrimes.View.CopyFromCPU(hostPrimes);
var devicePow2 = accelerator.Allocate1D<ulong>(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<uint>(hostLastOne.Length);
deviceLastOne.View.CopyFromCPU(hostLastOne);
var deviceLastSeven = accelerator.Allocate1D<uint>(hostLastSeven.Length);
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);
}
}
PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs
+3 -4
View File
@@ -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<uint, Stride1D.Dense> smallPrimesView = default;
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2View = default;
ResiduePrimeViews primeViews = default;
if (_kernelType == GpuKernelType.Pow2Mod)
{
(smallPrimesView, smallPrimesPow2View) = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
}
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);
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
{
PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs
+2 -2
View File
@@ -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);
PerfectNumbers.Core/PrimesGenerator.cs
+98 -65
View File
@@ -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<uint> 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;
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastSeven(uint prime)
{
uint mod10 = prime % 10U;
return mod10 == 3U || mod10 == 7U || mod10 == 9U || prime == 11U;
}
}