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

Merge pull request #61 from mlemanczyk/codex/analyze-and-optimize-evenperfectbitscanner

Browse Source 
Use Montgomery modular multiply in by-divisor GPU tester
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-22 13:22:24 +02:00
2 parents 8ed8e56ec0 a6571c2487
commit e7668874a6
3 files changed
+611 -160
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

EvenPerfectBitScanner.Benchmarks/GpuModularArithmeticBenchmarks.cs
+243
View File
@@ -0,0 +1,243 @@
using System.Numerics;
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser]
[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
public class GpuModularArithmeticBenchmarks
{
private const int MaxBatchSize = 512;
private readonly ulong[] _exponents = new ulong[MaxBatchSize];
private readonly ulong[] _moduli = new ulong[MaxBatchSize];
private readonly ulong[] _mersenneModuli = new ulong[MaxBatchSize];
private readonly int[] _mersenneBitWidths = new int[MaxBatchSize];
private readonly Random _random = new(7);
[Params(32, 256)]
public int BatchSize { get; set; }
[GlobalSetup]
public void Setup()
{
for (int i = 0; i < MaxBatchSize; i++)
{
_exponents[i] = NextExponent();
_moduli[i] = NextOddModulus();
int bits = _random.Next(8, 62);
_mersenneBitWidths[i] = bits;
_mersenneModuli[i] = (1UL << bits) - 1UL;
}
}
[Benchmark(Baseline = true)]
public ulong BaselineMulMod()
{
ulong checksum = 0UL;
for (int i = 0; i < BatchSize; i++)
{
checksum ^= Pow2ModBaseline(_exponents[i], _moduli[i]);
}
return checksum;
}
[Benchmark]
public ulong MultiplyHighMulMod()
{
ulong checksum = 0UL;
for (int i = 0; i < BatchSize; i++)
{
checksum ^= Pow2ModMultiplyHigh(_exponents[i], _moduli[i]);
}
return checksum;
}
[Benchmark]
public ulong MersenneFoldingMulMod()
{
ulong checksum = 0UL;
for (int i = 0; i < BatchSize; i++)
{
checksum ^= Pow2ModMersenne(_exponents[i], _mersenneBitWidths[i], _mersenneModuli[i]);
}
return checksum;
}
private ulong NextExponent()
{
return (ulong)_random.Next(24, 63);
}
private ulong NextOddModulus()
{
ulong value = (ulong)_random.NextInt64(5L, long.MaxValue);
if ((value & 1UL) == 0UL)
{
value++;
}
return value;
}
private static ulong Pow2ModBaseline(ulong exponent, ulong modulus)
{
if (modulus <= 1UL)
{
return 0UL;
}
ulong result = 1UL % modulus;
ulong baseVal = 2UL % modulus;
ulong exp = exponent;
while (exp > 0UL)
{
if ((exp & 1UL) != 0UL)
{
result = MulModBaseline(result, baseVal, modulus);
}
exp >>= 1;
if (exp == 0UL)
{
break;
}
baseVal = MulModBaseline(baseVal, baseVal, modulus);
}
return result;
}
private static ulong Pow2ModMultiplyHigh(ulong exponent, ulong modulus)
{
if (modulus <= 1UL)
{
return 0UL;
}
ulong result = 1UL % modulus;
ulong baseVal = 2UL % modulus;
ulong exp = exponent;
while (exp > 0UL)
{
if ((exp & 1UL) != 0UL)
{
result = MulModMultiplyHigh(result, baseVal, modulus);
}
exp >>= 1;
if (exp == 0UL)
{
break;
}
baseVal = MulModMultiplyHigh(baseVal, baseVal, modulus);
}
return result;
}
private static ulong Pow2ModMersenne(ulong exponent, int bitWidth, ulong modulus)
{
if (modulus <= 1UL)
{
return 0UL;
}
ulong result = 1UL & modulus;
ulong baseVal = 2UL & modulus;
ulong exp = exponent;
while (exp > 0UL)
{
if ((exp & 1UL) != 0UL)
{
result = MulModMersenne(result, baseVal, modulus, bitWidth);
}
exp >>= 1;
if (exp == 0UL)
{
break;
}
baseVal = MulModMersenne(baseVal, baseVal, modulus, bitWidth);
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong MulModBaseline(ulong a, ulong b, ulong modulus)
{
if (modulus == 0UL)
{
return 0UL;
}
ulong result = 0UL;
ulong x = a % modulus;
ulong y = b;
while (y > 0UL)
{
if ((y & 1UL) != 0UL)
{
result += x;
if (result >= modulus)
{
result -= modulus;
}
}
x <<= 1;
if (x >= modulus)
{
x -= modulus;
}
y >>= 1;
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong MulModMultiplyHigh(ulong a, ulong b, ulong modulus)
{
if (modulus == 0UL)
{
return 0UL;
}
UInt128 product = (UInt128)a * b;
return (ulong)(product % modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong MulModMersenne(ulong a, ulong b, ulong modulus, int bitWidth)
{
UInt128 product = (UInt128)a * b;
ulong mask = modulus;
ulong low = (ulong)product & mask;
ulong high = (ulong)(product >> bitWidth);
ulong result = low + high;
while (result > mask)
{
result = (result & mask) + (result >> bitWidth);
}
if (result >= modulus)
{
result -= modulus;
}
return result;
}
}
EvenPerfectBitScanner/Program.cs
+195 -53
View File
@@ -59,6 +59,7 @@ internal static class Program
private const int ByDivisorStateComposite = 1; private const int ByDivisorStateComposite = 1;
private const int ByDivisorStateCompleted = 2; private const int ByDivisorStateCompleted = 2;
private const int ByDivisorStateCompletedDetailed = 3; private const int ByDivisorStateCompletedDetailed = 3;
private const int DivisorAllocationBlockSize = 64;
private struct ByDivisorPrimeState private struct ByDivisorPrimeState
{ {
@@ -763,7 +764,7 @@ internal static class Program
return; return;
} }
List<ByDivisorPrimeState> states = new(primes.Count); List<ByDivisorPrimeState> states = new(primes.Count);
_byDivisorPrecheckOnly = true; _byDivisorPrecheckOnly = true;
foreach (ulong prime in primes) foreach (ulong prime in primes)
@@ -799,6 +800,8 @@ internal static class Program
return; return;
} }
states.Sort((left, right) => left.AllowedMax.CompareTo(right.AllowedMax));
int stateCount = states.Count; int stateCount = states.Count;
ulong[] primeValues = new ulong[stateCount]; ulong[] primeValues = new ulong[stateCount];
ulong[] allowedMaxValues = new ulong[stateCount]; ulong[] allowedMaxValues = new ulong[stateCount];
@@ -809,6 +812,7 @@ internal static class Program
var state = states[i]; var state = states[i];
primeValues[i] = state.Prime; primeValues[i] = state.Prime;
allowedMaxValues[i] = state.AllowedMax; allowedMaxValues[i] = state.AllowedMax;
stateFlags[i] = ByDivisorStateActive;
} }
states.Clear(); states.Clear();
@@ -820,6 +824,15 @@ internal static class Program
int finalizerState = 0; int finalizerState = 0;
int finalizationCompleted = 0; int finalizationCompleted = 0;
int remainingStates = stateCount; int remainingStates = stateCount;
int activeStartIndex = 0;
long[] activeStateMask = new long[(stateCount + 63) >> 6];
for (int i = 0; i < stateCount; i++)
{
int wordIndex = i >> 6;
int bitIndex = i & 63;
activeStateMask[wordIndex] |= 1L << bitIndex;
}
Task[] workers = new Task[Math.Max(1, threadCount)]; Task[] workers = new Task[Math.Max(1, threadCount)];
@@ -831,8 +844,12 @@ internal static class Program
byte[] hitsBuffer = ArrayPool<byte>.Shared.Rent(stateCount); byte[] hitsBuffer = ArrayPool<byte>.Shared.Rent(stateCount);
ulong[] primeBuffer = ArrayPool<ulong>.Shared.Rent(stateCount); ulong[] primeBuffer = ArrayPool<ulong>.Shared.Rent(stateCount);
int[] indexBuffer = ArrayPool<int>.Shared.Rent(stateCount); int[] indexBuffer = ArrayPool<int>.Shared.Rent(stateCount);
List<PendingResult> completions = new(stateCount); PendingResult[] completionsBuffer = ArrayPool<PendingResult>.Shared.Rent(stateCount);
List<PendingResult> composites = new(stateCount); PendingResult[] compositesBuffer = ArrayPool<PendingResult>.Shared.Rent(stateCount);
int completionsCount = 0;
int compositesCount = 0;
ulong localDivisorCursor = 0UL;
int localDivisorsRemaining = 0;
try try
{ {
@@ -847,7 +864,7 @@ internal static class Program
} }
bool exhausted; bool exhausted;
ulong divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted); ulong divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining);
if (divisor == 0UL) if (divisor == 0UL)
{ {
@@ -864,17 +881,11 @@ internal static class Program
if (Interlocked.CompareExchange(ref finalizerState, 1, 0) == 0) if (Interlocked.CompareExchange(ref finalizerState, 1, 0) == 0)
{ {
completions.Clear(); FinalizeRemainingStates(primeValues, allowedMaxValues, stateFlags, ref remainingStates, ref finalDivisorBits, completionsBuffer, ref completionsCount, activeStateMask);
FinalizeRemainingStates(primeValues, allowedMaxValues, stateFlags, ref remainingStates, ref finalDivisorBits, completions);
if (completions.Count > 0) if (completionsCount > 0)
{ {
foreach (PendingResult result in completions) FlushPendingResults(completionsBuffer, ref completionsCount);
{
PrintResult(result.Prime, searchedMersenne: true, detailedCheck: result.DetailedCheck, passedAllTests: result.PassedAllTests);
}
completions.Clear();
} }
Volatile.Write(ref finalizationCompleted, 1); Volatile.Write(ref finalizationCompleted, 1);
@@ -895,16 +906,11 @@ internal static class Program
continue; continue;
} }
int activeCount = BuildPrimeBuffer(divisor, primeValues, allowedMaxValues, stateFlags, primeBuffer, indexBuffer, completions, ref remainingStates); int activeCount = BuildPrimeBuffer(divisor, primeValues, allowedMaxValues, stateFlags, primeBuffer, indexBuffer, completionsBuffer, ref completionsCount, ref remainingStates, activeStateMask, ref activeStartIndex);
if (completions.Count > 0) if (completionsCount > 0)
{ {
foreach (PendingResult result in completions) FlushPendingResults(completionsBuffer, ref completionsCount);
{
PrintResult(result.Prime, searchedMersenne: true, detailedCheck: result.DetailedCheck, passedAllTests: result.PassedAllTests);
}
completions.Clear();
} }
if (activeCount == 0) if (activeCount == 0)
@@ -926,24 +932,37 @@ internal static class Program
int index = indexBuffer[j]; int index = indexBuffer[j];
if (Interlocked.CompareExchange(ref stateFlags[index], ByDivisorStateComposite, ByDivisorStateActive) == ByDivisorStateActive) if (Interlocked.CompareExchange(ref stateFlags[index], ByDivisorStateComposite, ByDivisorStateActive) == ByDivisorStateActive)
{ {
ClearActiveMask(activeStateMask, index);
Interlocked.Decrement(ref remainingStates); Interlocked.Decrement(ref remainingStates);
composites.Add(new PendingResult(primeValues[index], detailedCheck: true, passedAllTests: false)); compositesBuffer[compositesCount++] = new PendingResult(primeValues[index], detailedCheck: true, passedAllTests: false);
} }
}
if (composites.Count > 0) if (compositesCount == compositesBuffer.Length)
{
foreach (PendingResult result in composites)
{ {
PrintResult(result.Prime, searchedMersenne: true, detailedCheck: result.DetailedCheck, passedAllTests: result.PassedAllTests); FlushPendingResults(compositesBuffer, ref compositesCount);
} }
composites.Clear();
} }
if (compositesCount > 0)
{
FlushPendingResults(compositesBuffer, ref compositesCount);
}
}
if (completionsCount > 0)
{
FlushPendingResults(completionsBuffer, ref completionsCount);
}
if (compositesCount > 0)
{
FlushPendingResults(compositesBuffer, ref compositesCount);
} }
} }
finally finally
{ {
ArrayPool<PendingResult>.Shared.Return(compositesBuffer, clearArray: true);
ArrayPool<PendingResult>.Shared.Return(completionsBuffer, clearArray: true);
ArrayPool<int>.Shared.Return(indexBuffer, clearArray: true); ArrayPool<int>.Shared.Return(indexBuffer, clearArray: true);
ArrayPool<ulong>.Shared.Return(primeBuffer, clearArray: true); ArrayPool<ulong>.Shared.Return(primeBuffer, clearArray: true);
ArrayPool<byte>.Shared.Return(hitsBuffer, clearArray: true); ArrayPool<byte>.Shared.Return(hitsBuffer, clearArray: true);
@@ -954,10 +973,19 @@ internal static class Program
Task.WaitAll(workers); Task.WaitAll(workers);
} }
private static ulong AcquireNextDivisor(ref ulong nextDivisor, ulong divisorLimit, ref int divisorsExhaustedFlag, ref long finalDivisorBits, out bool exhausted) private static ulong AcquireNextDivisor(ref ulong nextDivisor, ulong divisorLimit, ref int divisorsExhaustedFlag, ref long finalDivisorBits, out bool exhausted, ref ulong localDivisorCursor, ref int localDivisorsRemaining)
{ {
ref long nextDivisorBits = ref Unsafe.As<ulong, long>(ref nextDivisor); ref long nextDivisorBits = ref Unsafe.As<ulong, long>(ref nextDivisor);
if (localDivisorsRemaining > 0)
{
ulong currentValue = localDivisorCursor;
localDivisorCursor += 2UL;
localDivisorsRemaining--;
exhausted = false;
return currentValue;
}
while (true) while (true)
{ {
if (Volatile.Read(ref divisorsExhaustedFlag) != 0) if (Volatile.Read(ref divisorsExhaustedFlag) != 0)
@@ -979,47 +1007,124 @@ internal static class Program
return 0UL; return 0UL;
} }
ulong nextValue = currentValue + 2UL; ulong maximumNext = divisorLimit >= ulong.MaxValue - 1UL ? ulong.MaxValue : divisorLimit + 2UL;
long nextBits = unchecked((long)nextValue); ulong blockStride = unchecked((ulong)(DivisorAllocationBlockSize * 2));
if (Interlocked.CompareExchange(ref nextDivisorBits, nextBits, currentBits) == currentBits) ulong requestedNext = currentValue > ulong.MaxValue - blockStride ? ulong.MaxValue : currentValue + blockStride;
if (requestedNext > maximumNext)
{ {
exhausted = false; requestedNext = maximumNext;
return currentValue;
} }
long nextBits = unchecked((long)requestedNext);
if (Interlocked.CompareExchange(ref nextDivisorBits, nextBits, currentBits) != currentBits)
{
continue;
}
ulong available = requestedNext - currentValue;
if (available == 0UL)
{
continue;
}
int count = (int)(available >> 1);
if (count <= 0)
{
continue;
}
localDivisorCursor = currentValue + 2UL;
localDivisorsRemaining = count - 1;
exhausted = false;
return currentValue;
} }
} }
private static int BuildPrimeBuffer(ulong divisor, ulong[] primeValues, ulong[] allowedMaxValues, int[] stateFlags, ulong[] primeBuffer, int[] indexBuffer, List<PendingResult> completions, ref int remainingStates) private static int BuildPrimeBuffer(ulong divisor, ulong[] primeValues, ulong[] allowedMaxValues, int[] stateFlags, ulong[] primeBuffer, int[] indexBuffer, PendingResult[] completionsBuffer, ref int completionsCount, ref int remainingStates, long[] activeStateMask, ref int activeStartIndex)
{ {
int activeCount = 0; int length = primeValues.Length;
int startIndex = Volatile.Read(ref activeStartIndex);
int index = startIndex;
for (int i = 0; i < primeValues.Length; i++) while (index < length && allowedMaxValues[index] < divisor)
{ {
if (allowedMaxValues[i] < divisor) if (Interlocked.CompareExchange(ref stateFlags[index], ByDivisorStateCompletedDetailed, ByDivisorStateActive) == ByDivisorStateActive)
{ {
if (Interlocked.CompareExchange(ref stateFlags[i], ByDivisorStateCompletedDetailed, ByDivisorStateActive) == ByDivisorStateActive) ClearActiveMask(activeStateMask, index);
Interlocked.Decrement(ref remainingStates);
if (completionsCount == completionsBuffer.Length)
{ {
Interlocked.Decrement(ref remainingStates); FlushPendingResults(completionsBuffer, ref completionsCount);
completions.Add(new PendingResult(primeValues[i], detailedCheck: true, passedAllTests: true)); }
completionsBuffer[completionsCount++] = new PendingResult(primeValues[index], detailedCheck: true, passedAllTests: true);
}
else
{
ClearActiveMask(activeStateMask, index);
}
index++;
}
if (index != startIndex)
{
Volatile.Write(ref activeStartIndex, index);
}
int activeCount = 0;
while (index < length)
{
int wordIndex = index >> 6;
ulong word = unchecked((ulong)Volatile.Read(ref activeStateMask[wordIndex]));
if (word == 0UL)
{
index = (wordIndex + 1) << 6;
continue;
}
int bitOffset = index & 63;
if (bitOffset != 0)
{
ulong mask = ulong.MaxValue << bitOffset;
word &= mask;
if (word == 0UL)
{
index = (wordIndex + 1) << 6;
continue;
}
}
while (word != 0UL)
{
int bit = BitOperations.TrailingZeroCount(word);
int candidateIndex = (wordIndex << 6) + bit;
if (candidateIndex >= length)
{
return activeCount;
} }
continue; int state = Volatile.Read(ref stateFlags[candidateIndex]);
if (state == ByDivisorStateActive)
{
primeBuffer[activeCount] = primeValues[candidateIndex];
indexBuffer[activeCount] = candidateIndex;
activeCount++;
}
else
{
ClearActiveMask(activeStateMask, candidateIndex);
}
word &= word - 1UL;
} }
if (Volatile.Read(ref stateFlags[i]) != ByDivisorStateActive) index = (wordIndex + 1) << 6;
{
continue;
}
primeBuffer[activeCount] = primeValues[i];
indexBuffer[activeCount] = i;
activeCount++;
} }
return activeCount; return activeCount;
} }
private static void FinalizeRemainingStates(ulong[] primeValues, ulong[] allowedMaxValues, int[] stateFlags, ref int remainingStates, ref long finalDivisorBits, List<PendingResult> completions) private static void FinalizeRemainingStates(ulong[] primeValues, ulong[] allowedMaxValues, int[] stateFlags, ref int remainingStates, ref long finalDivisorBits, PendingResult[] completionsBuffer, ref int completionsCount, long[] activeStateMask)
{ {
ulong finalDivisorValue = unchecked((ulong)Volatile.Read(ref finalDivisorBits)); ulong finalDivisorValue = unchecked((ulong)Volatile.Read(ref finalDivisorBits));
@@ -1038,8 +1143,45 @@ internal static class Program
continue; continue;
} }
ClearActiveMask(activeStateMask, i);
Interlocked.Decrement(ref remainingStates); Interlocked.Decrement(ref remainingStates);
completions.Add(new PendingResult(primeValues[i], detailedCheck: detailed, passedAllTests: true)); if (completionsCount == completionsBuffer.Length)
{
FlushPendingResults(completionsBuffer, ref completionsCount);
}
completionsBuffer[completionsCount++] = new PendingResult(primeValues[i], detailedCheck: detailed, passedAllTests: true);
}
}
private static void FlushPendingResults(PendingResult[] buffer, ref int count)
{
for (int i = 0; i < count; i++)
{
PendingResult result = buffer[i];
PrintResult(result.Prime, searchedMersenne: true, detailedCheck: result.DetailedCheck, passedAllTests: result.PassedAllTests);
}
count = 0;
}
private static void ClearActiveMask(long[] activeStateMask, int index)
{
int wordIndex = index >> 6;
long bit = 1L << (index & 63);
while (true)
{
long current = Volatile.Read(ref activeStateMask[wordIndex]);
if ((current & bit) == 0)
{
return;
}
long next = current & ~bit;
if (Interlocked.CompareExchange(ref activeStateMask[wordIndex], next, current) == current)
{
return;
}
} }
} }
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+173 -107
View File
@@ -1,9 +1,10 @@
using System; using System;
using System.Buffers; using System.Buffers;
using System.Collections.Concurrent; using System.Collections.Concurrent;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices; using System.Runtime.InteropServices;
using ILGPU; using ILGPU;
using ILGPU.Algorithms;
using ILGPU.Runtime; using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Gpu;
@@ -16,13 +17,29 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
private bool _isConfigured; private bool _isConfigured;
private ulong _lastStatusDivisor; private ulong _lastStatusDivisor;
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelCache = new(); public readonly struct MontgomeryDivisorData
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelByPrimeCache = new(); {
public readonly ulong Modulus;
public readonly ulong NPrime;
public readonly ulong MontgomeryOne;
public readonly ulong MontgomeryTwo;
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetKernel(Accelerator accelerator) => public MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo)
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernel)); {
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetKernelByPrime(Accelerator accelerator) => Modulus = modulus;
_kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernelByPrime)); NPrime = nPrime;
MontgomeryOne = montgomeryOne;
MontgomeryTwo = montgomeryTwo;
}
}
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>> _kernelCache = new();
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<byte>>> _kernelByPrimeCache = new();
private Action<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>(CheckKernel));
private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<byte>> GetKernelByPrime(Accelerator accelerator) =>
_kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<byte>>(CheckKernelByPrime));
public int GpuBatchSize public int GpuBatchSize
{ {
@@ -90,15 +107,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
var kernel = GetKernel(accelerator); var kernel = GetKernel(accelerator);
int batchCapacity = _gpuBatchSize; int batchCapacity = _gpuBatchSize;
var divisorsBuffer = accelerator.Allocate1D<ulong>(batchCapacity); var divisorsBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(batchCapacity);
var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity); var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity);
ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity); ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity); byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
MontgomeryDivisorData[] divisorData = ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity);
bool composite = false; bool composite = false;
bool processedAll = false; bool processedAll = false;
ulong divisor = 3UL; ulong divisor = 3UL;
try try
{ {
@@ -135,17 +153,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
break; break;
} }
var divisorView = divisorsBuffer.View.SubView(0, batchSize); Span<MontgomeryDivisorData> divisorDataSpan = divisorData.AsSpan(0, batchSize);
var hitsView = hitsBuffer.View.SubView(0, batchSize); var divisorView = divisorsBuffer.View.SubView(0, batchSize);
var hitsView = hitsBuffer.View.SubView(0, batchSize);
Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize); Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize);
Span<byte> hitsSpan = hits.AsSpan(0, batchSize); Span<byte> hitsSpan = hits.AsSpan(0, batchSize);
divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize); for (int i = 0; i < batchSize; i++)
kernel(batchSize, prime, divisorView, hitsView); {
accelerator.Synchronize(); divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]);
}
hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize);
kernel(batchSize, prime, divisorView, hitsView);
accelerator.Synchronize();
hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize);
for (int i = 0; i < batchSize; i++) for (int i = 0; i < batchSize; i++)
{ {
@@ -174,12 +198,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
} }
finally finally
{ {
ArrayPool<byte>.Shared.Return(hits, clearArray: true); ArrayPool<byte>.Shared.Return(hits, clearArray: true);
ArrayPool<ulong>.Shared.Return(divisors, clearArray: true); ArrayPool<ulong>.Shared.Return(divisors, clearArray: true);
hitsBuffer.Dispose(); ArrayPool<MontgomeryDivisorData>.Shared.Return(divisorData, clearArray: true);
divisorsBuffer.Dispose(); hitsBuffer.Dispose();
gpuLease.Dispose(); divisorsBuffer.Dispose();
} gpuLease.Dispose();
}
return composite; return composite;
} }
@@ -224,25 +249,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
return Math.Min(maxByPrime, _divisorLimit); return Math.Min(maxByPrime, _divisorLimit);
} }
private static void CheckKernel(Index1D index, ulong prime, ArrayView<ulong> divisors, ArrayView<byte> hits) private static void CheckKernel(Index1D index, ulong prime, ArrayView<MontgomeryDivisorData> divisors, ArrayView<byte> hits)
{ {
ulong divisor = divisors[index]; MontgomeryDivisorData divisor = divisors[index];
if (divisor <= 1UL || (divisor & 1UL) == 0UL) if (divisor.Modulus <= 1UL || (divisor.Modulus & 1UL) == 0UL)
{ {
hits[index] = 0; hits[index] = 0;
return; return;
} }
ulong pow = Pow2Mod(prime, divisor); ulong pow = Pow2Mod(prime, divisor);
hits[index] = pow == 1UL ? (byte)1 : (byte)0; hits[index] = pow == 1UL ? (byte)1 : (byte)0;
} }
public sealed class DivisorScanSession : IDisposable public sealed class DivisorScanSession : IDisposable
{ {
private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
private readonly GpuContextPool.GpuContextLease _lease; private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator; private readonly Accelerator _accelerator;
private readonly Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> _kernel; private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<byte>> _kernel;
private readonly MemoryBuffer1D<ulong, Stride1D.Dense> _primesBuffer; private readonly MemoryBuffer1D<ulong, Stride1D.Dense> _primesBuffer;
private readonly MemoryBuffer1D<byte, Stride1D.Dense> _hitsBuffer; private readonly MemoryBuffer1D<byte, Stride1D.Dense> _hitsBuffer;
private readonly ulong[] _primesHost; private readonly ulong[] _primesHost;
@@ -254,7 +279,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
_owner = owner; _owner = owner;
_lease = GpuContextPool.RentPreferred(preferCpu: false); _lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator; _accelerator = _lease.Accelerator;
_kernel = owner.GetKernelByPrime(_accelerator); _kernel = owner.GetKernelByPrime(_accelerator);
_primesBuffer = _accelerator.Allocate1D<ulong>(owner._gpuBatchSize); _primesBuffer = _accelerator.Allocate1D<ulong>(owner._gpuBatchSize);
_hitsBuffer = _accelerator.Allocate1D<byte>(owner._gpuBatchSize); _hitsBuffer = _accelerator.Allocate1D<byte>(owner._gpuBatchSize);
_primesHost = ArrayPool<ulong>.Shared.Rent(owner._gpuBatchSize); _primesHost = ArrayPool<ulong>.Shared.Rent(owner._gpuBatchSize);
@@ -275,8 +300,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
int offset = 0; int offset = 0;
int consoleStatus = 0; int consoleStatus = 0;
while (offset < primes.Length) MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
{ while (offset < primes.Length)
{
int batchSize = Math.Min(_owner._gpuBatchSize, primes.Length - offset); int batchSize = Math.Min(_owner._gpuBatchSize, primes.Length - offset);
Span<ulong> primeSpan = _primesHost.AsSpan(0, batchSize); Span<ulong> primeSpan = _primesHost.AsSpan(0, batchSize);
primes.Slice(offset, batchSize).CopyTo(primeSpan); primes.Slice(offset, batchSize).CopyTo(primeSpan);
@@ -291,7 +317,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
var hitsView = _hitsBuffer.View.SubView(0, batchSize); var hitsView = _hitsBuffer.View.SubView(0, batchSize);
primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primeSpan), batchSize); primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primeSpan), batchSize);
_kernel(batchSize, divisor, primesView, hitsView); _kernel(batchSize, divisorData, primesView, hitsView);
_accelerator.Synchronize(); _accelerator.Synchronize();
Span<byte> hitSpan = _hitsHost.AsSpan(0, batchSize); Span<byte> hitSpan = _hitsHost.AsSpan(0, batchSize);
@@ -317,80 +343,120 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
} }
} }
private static ulong Pow2Mod(ulong exponent, ulong modulus) private static ulong Pow2Mod(ulong exponent, MontgomeryDivisorData divisor)
{ {
if (modulus <= 1UL) ulong modulus = divisor.Modulus;
{ if (modulus <= 1UL || (modulus & 1UL) == 0UL)
return 0UL; {
} return 0UL;
}
ulong result = 1UL % modulus; ulong result = divisor.MontgomeryOne;
ulong baseVal = 2UL % modulus; ulong baseVal = divisor.MontgomeryTwo;
ulong exp = exponent; ulong exp = exponent;
while (exp > 0UL) while (exp > 0UL)
{ {
if ((exp & 1UL) != 0UL) if ((exp & 1UL) != 0UL)
{ {
result = MulMod(result, baseVal, modulus); result = MontgomeryMultiply(result, baseVal, modulus, divisor.NPrime);
} }
exp >>= 1; exp >>= 1;
if (exp == 0UL) if (exp == 0UL)
{ {
break; break;
} }
baseVal = MulMod(baseVal, baseVal, modulus); baseVal = MontgomeryMultiply(baseVal, baseVal, modulus, divisor.NPrime);
} }
return result; return MontgomeryMultiply(result, 1UL, modulus, divisor.NPrime);
} }
private static ulong MulMod(ulong a, ulong b, ulong modulus) [MethodImpl(MethodImplOptions.AggressiveInlining)]
{ private static ulong MontgomeryMultiply(ulong a, ulong b, ulong modulus, ulong nPrime)
if (modulus == 0UL) {
{ ulong tLow = unchecked(a * b);
return 0UL; ulong tHigh = MultiplyHigh(a, b);
} ulong m = unchecked(tLow * nPrime);
ulong mnLow = unchecked(m * modulus);
ulong mnHigh = MultiplyHigh(m, modulus);
ulong result = 0UL; ulong sumLow = unchecked(tLow + mnLow);
ulong x = a % modulus; ulong carry = sumLow < tLow ? 1UL : 0UL;
ulong y = b; ulong sumHigh = unchecked(tHigh + mnHigh + carry);
while (y > 0UL)
{
if ((y & 1UL) != 0UL)
{
result += x;
if (result >= modulus)
{
result -= modulus;
}
}
x <<= 1; ulong result = sumHigh;
if (x >= modulus) if (result >= modulus)
{ {
x -= modulus; result -= modulus;
} }
y >>= 1; return result;
} }
return result; private static void CheckKernelByPrime(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> primes, ArrayView<byte> hits)
} {
ulong prime = primes[index];
if (divisor.Modulus <= 1UL || (divisor.Modulus & 1UL) == 0UL)
{
hits[index] = 0;
return;
}
private static void CheckKernelByPrime(Index1D index, ulong divisor, ArrayView<ulong> primes, ArrayView<byte> hits) ulong pow = Pow2Mod(prime, divisor);
{ hits[index] = pow == 1UL ? (byte)1 : (byte)0;
ulong prime = primes[index]; }
if (divisor <= 1UL || (divisor & 1UL) == 0UL)
{
hits[index] = 0;
return;
}
ulong pow = Pow2Mod(prime, divisor); private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
hits[index] = pow == 1UL ? (byte)1 : (byte)0; {
} if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
}
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
return new MontgomeryDivisorData(modulus, nPrime, montgomeryOne, montgomeryTwo);
}
private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus)
{
UInt128 r = UInt128.One << 64;
UInt128 result = ((UInt128)value * r) % modulus;
return (ulong)result;
}
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong MultiplyHigh(ulong x, ulong y)
{
ulong xLow = (uint)x;
ulong xHigh = x >> 32;
ulong yLow = (uint)y;
ulong yHigh = y >> 32;
ulong lowProduct = xLow * yLow;
ulong mid = xHigh * yLow + (lowProduct >> 32);
ulong midLow = mid & 0xFFFFFFFFUL;
ulong midHigh = mid >> 32;
midLow += xLow * yHigh;
ulong resultHigh = xHigh * yHigh + midHigh + (midLow >> 32);
return resultHigh;
}
public ulong DivisorLimit public ulong DivisorLimit
{ {