Accelerate by-divisor GPU tester with Montgomery math
This commit is contained in:
Marek Lemańczyk committed
1 parent
8ed8e56ec0
commit
a6571c2487
3 files changed
+611
-160
No files matched your search
@@ -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;
|
||||
}
|
||||
}
|
||||
@@ -59,6 +59,7 @@ internal static class Program
|
||||
private const int ByDivisorStateComposite = 1;
|
||||
private const int ByDivisorStateCompleted = 2;
|
||||
private const int ByDivisorStateCompletedDetailed = 3;
|
||||
private const int DivisorAllocationBlockSize = 64;
|
||||
|
||||
private struct ByDivisorPrimeState
|
||||
{
|
||||
@@ -763,7 +764,7 @@ internal static class Program
|
||||
return;
|
||||
}
|
||||
|
||||
List<ByDivisorPrimeState> states = new(primes.Count);
|
||||
List<ByDivisorPrimeState> states = new(primes.Count);
|
||||
|
||||
_byDivisorPrecheckOnly = true;
|
||||
foreach (ulong prime in primes)
|
||||
@@ -799,6 +800,8 @@ internal static class Program
|
||||
return;
|
||||
}
|
||||
|
||||
states.Sort((left, right) => left.AllowedMax.CompareTo(right.AllowedMax));
|
||||
|
||||
int stateCount = states.Count;
|
||||
ulong[] primeValues = new ulong[stateCount];
|
||||
ulong[] allowedMaxValues = new ulong[stateCount];
|
||||
@@ -809,6 +812,7 @@ internal static class Program
|
||||
var state = states[i];
|
||||
primeValues[i] = state.Prime;
|
||||
allowedMaxValues[i] = state.AllowedMax;
|
||||
stateFlags[i] = ByDivisorStateActive;
|
||||
}
|
||||
|
||||
states.Clear();
|
||||
@@ -820,6 +824,15 @@ internal static class Program
|
||||
int finalizerState = 0;
|
||||
int finalizationCompleted = 0;
|
||||
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)];
|
||||
|
||||
@@ -831,8 +844,12 @@ internal static class Program
|
||||
byte[] hitsBuffer = ArrayPool<byte>.Shared.Rent(stateCount);
|
||||
ulong[] primeBuffer = ArrayPool<ulong>.Shared.Rent(stateCount);
|
||||
int[] indexBuffer = ArrayPool<int>.Shared.Rent(stateCount);
|
||||
List<PendingResult> completions = new(stateCount);
|
||||
List<PendingResult> composites = new(stateCount);
|
||||
PendingResult[] completionsBuffer = ArrayPool<PendingResult>.Shared.Rent(stateCount);
|
||||
PendingResult[] compositesBuffer = ArrayPool<PendingResult>.Shared.Rent(stateCount);
|
||||
int completionsCount = 0;
|
||||
int compositesCount = 0;
|
||||
ulong localDivisorCursor = 0UL;
|
||||
int localDivisorsRemaining = 0;
|
||||
|
||||
try
|
||||
{
|
||||
@@ -847,7 +864,7 @@ internal static class Program
|
||||
}
|
||||
|
||||
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)
|
||||
{
|
||||
@@ -864,17 +881,11 @@ internal static class Program
|
||||
|
||||
if (Interlocked.CompareExchange(ref finalizerState, 1, 0) == 0)
|
||||
{
|
||||
completions.Clear();
|
||||
FinalizeRemainingStates(primeValues, allowedMaxValues, stateFlags, ref remainingStates, ref finalDivisorBits, completions);
|
||||
FinalizeRemainingStates(primeValues, allowedMaxValues, stateFlags, ref remainingStates, ref finalDivisorBits, completionsBuffer, ref completionsCount, activeStateMask);
|
||||
|
||||
if (completions.Count > 0)
|
||||
if (completionsCount > 0)
|
||||
{
|
||||
foreach (PendingResult result in completions)
|
||||
{
|
||||
PrintResult(result.Prime, searchedMersenne: true, detailedCheck: result.DetailedCheck, passedAllTests: result.PassedAllTests);
|
||||
}
|
||||
|
||||
completions.Clear();
|
||||
FlushPendingResults(completionsBuffer, ref completionsCount);
|
||||
}
|
||||
|
||||
Volatile.Write(ref finalizationCompleted, 1);
|
||||
@@ -895,16 +906,11 @@ internal static class Program
|
||||
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)
|
||||
{
|
||||
PrintResult(result.Prime, searchedMersenne: true, detailedCheck: result.DetailedCheck, passedAllTests: result.PassedAllTests);
|
||||
}
|
||||
|
||||
completions.Clear();
|
||||
FlushPendingResults(completionsBuffer, ref completionsCount);
|
||||
}
|
||||
|
||||
if (activeCount == 0)
|
||||
@@ -926,24 +932,37 @@ internal static class Program
|
||||
int index = indexBuffer[j];
|
||||
if (Interlocked.CompareExchange(ref stateFlags[index], ByDivisorStateComposite, ByDivisorStateActive) == ByDivisorStateActive)
|
||||
{
|
||||
ClearActiveMask(activeStateMask, index);
|
||||
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)
|
||||
{
|
||||
foreach (PendingResult result in composites)
|
||||
if (compositesCount == compositesBuffer.Length)
|
||||
{
|
||||
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
|
||||
{
|
||||
ArrayPool<PendingResult>.Shared.Return(compositesBuffer, clearArray: true);
|
||||
ArrayPool<PendingResult>.Shared.Return(completionsBuffer, clearArray: true);
|
||||
ArrayPool<int>.Shared.Return(indexBuffer, clearArray: true);
|
||||
ArrayPool<ulong>.Shared.Return(primeBuffer, clearArray: true);
|
||||
ArrayPool<byte>.Shared.Return(hitsBuffer, clearArray: true);
|
||||
@@ -954,10 +973,19 @@ internal static class Program
|
||||
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);
|
||||
|
||||
if (localDivisorsRemaining > 0)
|
||||
{
|
||||
ulong currentValue = localDivisorCursor;
|
||||
localDivisorCursor += 2UL;
|
||||
localDivisorsRemaining--;
|
||||
exhausted = false;
|
||||
return currentValue;
|
||||
}
|
||||
|
||||
while (true)
|
||||
{
|
||||
if (Volatile.Read(ref divisorsExhaustedFlag) != 0)
|
||||
@@ -979,47 +1007,124 @@ internal static class Program
|
||||
return 0UL;
|
||||
}
|
||||
|
||||
ulong nextValue = currentValue + 2UL;
|
||||
long nextBits = unchecked((long)nextValue);
|
||||
if (Interlocked.CompareExchange(ref nextDivisorBits, nextBits, currentBits) == currentBits)
|
||||
ulong maximumNext = divisorLimit >= ulong.MaxValue - 1UL ? ulong.MaxValue : divisorLimit + 2UL;
|
||||
ulong blockStride = unchecked((ulong)(DivisorAllocationBlockSize * 2));
|
||||
ulong requestedNext = currentValue > ulong.MaxValue - blockStride ? ulong.MaxValue : currentValue + blockStride;
|
||||
if (requestedNext > maximumNext)
|
||||
{
|
||||
exhausted = false;
|
||||
return currentValue;
|
||||
requestedNext = maximumNext;
|
||||
}
|
||||
|
||||
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);
|
||||
completions.Add(new PendingResult(primeValues[i], detailedCheck: true, passedAllTests: true));
|
||||
FlushPendingResults(completionsBuffer, ref completionsCount);
|
||||
}
|
||||
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)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
primeBuffer[activeCount] = primeValues[i];
|
||||
indexBuffer[activeCount] = i;
|
||||
activeCount++;
|
||||
index = (wordIndex + 1) << 6;
|
||||
}
|
||||
|
||||
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));
|
||||
|
||||
@@ -1038,8 +1143,45 @@ internal static class Program
|
||||
continue;
|
||||
}
|
||||
|
||||
ClearActiveMask(activeStateMask, i);
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -1,9 +1,10 @@
|
||||
using System;
|
||||
using System.Buffers;
|
||||
using System.Collections.Concurrent;
|
||||
using System.Numerics;
|
||||
using System.Runtime.CompilerServices;
|
||||
using System.Runtime.InteropServices;
|
||||
using ILGPU;
|
||||
using ILGPU.Algorithms;
|
||||
using ILGPU.Runtime;
|
||||
|
||||
namespace PerfectNumbers.Core.Gpu;
|
||||
@@ -16,13 +17,29 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
private bool _isConfigured;
|
||||
private ulong _lastStatusDivisor;
|
||||
|
||||
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelCache = new();
|
||||
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelByPrimeCache = new();
|
||||
public readonly struct MontgomeryDivisorData
|
||||
{
|
||||
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) =>
|
||||
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
|
||||
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetKernelByPrime(Accelerator accelerator) =>
|
||||
_kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernelByPrime));
|
||||
public MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo)
|
||||
{
|
||||
Modulus = modulus;
|
||||
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
|
||||
{
|
||||
@@ -90,15 +107,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
var kernel = GetKernel(accelerator);
|
||||
|
||||
int batchCapacity = _gpuBatchSize;
|
||||
var divisorsBuffer = accelerator.Allocate1D<ulong>(batchCapacity);
|
||||
var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity);
|
||||
var divisorsBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(batchCapacity);
|
||||
var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity);
|
||||
|
||||
ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
|
||||
byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
|
||||
ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
|
||||
byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
|
||||
MontgomeryDivisorData[] divisorData = ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity);
|
||||
|
||||
bool composite = false;
|
||||
bool processedAll = false;
|
||||
ulong divisor = 3UL;
|
||||
bool composite = false;
|
||||
bool processedAll = false;
|
||||
ulong divisor = 3UL;
|
||||
|
||||
try
|
||||
{
|
||||
@@ -135,17 +153,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
break;
|
||||
}
|
||||
|
||||
var divisorView = divisorsBuffer.View.SubView(0, batchSize);
|
||||
var hitsView = hitsBuffer.View.SubView(0, batchSize);
|
||||
Span<MontgomeryDivisorData> divisorDataSpan = divisorData.AsSpan(0, batchSize);
|
||||
var divisorView = divisorsBuffer.View.SubView(0, batchSize);
|
||||
var hitsView = hitsBuffer.View.SubView(0, batchSize);
|
||||
|
||||
Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize);
|
||||
Span<byte> hitsSpan = hits.AsSpan(0, batchSize);
|
||||
Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize);
|
||||
Span<byte> hitsSpan = hits.AsSpan(0, batchSize);
|
||||
|
||||
divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize);
|
||||
kernel(batchSize, prime, divisorView, hitsView);
|
||||
accelerator.Synchronize();
|
||||
for (int i = 0; i < batchSize; i++)
|
||||
{
|
||||
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++)
|
||||
{
|
||||
@@ -174,12 +198,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
}
|
||||
finally
|
||||
{
|
||||
ArrayPool<byte>.Shared.Return(hits, clearArray: true);
|
||||
ArrayPool<ulong>.Shared.Return(divisors, clearArray: true);
|
||||
hitsBuffer.Dispose();
|
||||
divisorsBuffer.Dispose();
|
||||
gpuLease.Dispose();
|
||||
}
|
||||
ArrayPool<byte>.Shared.Return(hits, clearArray: true);
|
||||
ArrayPool<ulong>.Shared.Return(divisors, clearArray: true);
|
||||
ArrayPool<MontgomeryDivisorData>.Shared.Return(divisorData, clearArray: true);
|
||||
hitsBuffer.Dispose();
|
||||
divisorsBuffer.Dispose();
|
||||
gpuLease.Dispose();
|
||||
}
|
||||
|
||||
return composite;
|
||||
}
|
||||
@@ -224,25 +249,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
return Math.Min(maxByPrime, _divisorLimit);
|
||||
}
|
||||
|
||||
private static void CheckKernel(Index1D index, ulong prime, ArrayView<ulong> divisors, ArrayView<byte> hits)
|
||||
{
|
||||
ulong divisor = divisors[index];
|
||||
if (divisor <= 1UL || (divisor & 1UL) == 0UL)
|
||||
{
|
||||
hits[index] = 0;
|
||||
return;
|
||||
}
|
||||
private static void CheckKernel(Index1D index, ulong prime, ArrayView<MontgomeryDivisorData> divisors, ArrayView<byte> hits)
|
||||
{
|
||||
MontgomeryDivisorData divisor = divisors[index];
|
||||
if (divisor.Modulus <= 1UL || (divisor.Modulus & 1UL) == 0UL)
|
||||
{
|
||||
hits[index] = 0;
|
||||
return;
|
||||
}
|
||||
|
||||
ulong pow = Pow2Mod(prime, divisor);
|
||||
hits[index] = pow == 1UL ? (byte)1 : (byte)0;
|
||||
}
|
||||
ulong pow = Pow2Mod(prime, divisor);
|
||||
hits[index] = pow == 1UL ? (byte)1 : (byte)0;
|
||||
}
|
||||
|
||||
public sealed class DivisorScanSession : IDisposable
|
||||
{
|
||||
private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
|
||||
private readonly GpuContextPool.GpuContextLease _lease;
|
||||
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<byte, Stride1D.Dense> _hitsBuffer;
|
||||
private readonly ulong[] _primesHost;
|
||||
@@ -254,7 +279,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
_owner = owner;
|
||||
_lease = GpuContextPool.RentPreferred(preferCpu: false);
|
||||
_accelerator = _lease.Accelerator;
|
||||
_kernel = owner.GetKernelByPrime(_accelerator);
|
||||
_kernel = owner.GetKernelByPrime(_accelerator);
|
||||
_primesBuffer = _accelerator.Allocate1D<ulong>(owner._gpuBatchSize);
|
||||
_hitsBuffer = _accelerator.Allocate1D<byte>(owner._gpuBatchSize);
|
||||
_primesHost = ArrayPool<ulong>.Shared.Rent(owner._gpuBatchSize);
|
||||
@@ -275,8 +300,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
|
||||
int offset = 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);
|
||||
Span<ulong> primeSpan = _primesHost.AsSpan(0, batchSize);
|
||||
primes.Slice(offset, batchSize).CopyTo(primeSpan);
|
||||
@@ -291,7 +317,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
var hitsView = _hitsBuffer.View.SubView(0, batchSize);
|
||||
|
||||
primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primeSpan), batchSize);
|
||||
_kernel(batchSize, divisor, primesView, hitsView);
|
||||
_kernel(batchSize, divisorData, primesView, hitsView);
|
||||
_accelerator.Synchronize();
|
||||
|
||||
Span<byte> hitSpan = _hitsHost.AsSpan(0, batchSize);
|
||||
@@ -317,80 +343,120 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
|
||||
}
|
||||
}
|
||||
|
||||
private static ulong Pow2Mod(ulong exponent, ulong modulus)
|
||||
{
|
||||
if (modulus <= 1UL)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
private static ulong Pow2Mod(ulong exponent, MontgomeryDivisorData divisor)
|
||||
{
|
||||
ulong modulus = divisor.Modulus;
|
||||
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
|
||||
ulong result = 1UL % modulus;
|
||||
ulong baseVal = 2UL % modulus;
|
||||
ulong exp = exponent;
|
||||
while (exp > 0UL)
|
||||
{
|
||||
if ((exp & 1UL) != 0UL)
|
||||
{
|
||||
result = MulMod(result, baseVal, modulus);
|
||||
}
|
||||
ulong result = divisor.MontgomeryOne;
|
||||
ulong baseVal = divisor.MontgomeryTwo;
|
||||
ulong exp = exponent;
|
||||
while (exp > 0UL)
|
||||
{
|
||||
if ((exp & 1UL) != 0UL)
|
||||
{
|
||||
result = MontgomeryMultiply(result, baseVal, modulus, divisor.NPrime);
|
||||
}
|
||||
|
||||
exp >>= 1;
|
||||
if (exp == 0UL)
|
||||
{
|
||||
break;
|
||||
}
|
||||
exp >>= 1;
|
||||
if (exp == 0UL)
|
||||
{
|
||||
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)
|
||||
{
|
||||
if (modulus == 0UL)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static ulong MontgomeryMultiply(ulong a, ulong b, ulong modulus, ulong nPrime)
|
||||
{
|
||||
ulong tLow = unchecked(a * b);
|
||||
ulong tHigh = MultiplyHigh(a, b);
|
||||
ulong m = unchecked(tLow * nPrime);
|
||||
ulong mnLow = unchecked(m * modulus);
|
||||
ulong mnHigh = MultiplyHigh(m, modulus);
|
||||
|
||||
ulong result = 0UL;
|
||||
ulong x = a % modulus;
|
||||
ulong y = b;
|
||||
while (y > 0UL)
|
||||
{
|
||||
if ((y & 1UL) != 0UL)
|
||||
{
|
||||
result += x;
|
||||
if (result >= modulus)
|
||||
{
|
||||
result -= modulus;
|
||||
}
|
||||
}
|
||||
ulong sumLow = unchecked(tLow + mnLow);
|
||||
ulong carry = sumLow < tLow ? 1UL : 0UL;
|
||||
ulong sumHigh = unchecked(tHigh + mnHigh + carry);
|
||||
|
||||
x <<= 1;
|
||||
if (x >= modulus)
|
||||
{
|
||||
x -= modulus;
|
||||
}
|
||||
ulong result = sumHigh;
|
||||
if (result >= 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 prime = primes[index];
|
||||
if (divisor <= 1UL || (divisor & 1UL) == 0UL)
|
||||
{
|
||||
hits[index] = 0;
|
||||
return;
|
||||
}
|
||||
ulong pow = Pow2Mod(prime, divisor);
|
||||
hits[index] = pow == 1UL ? (byte)1 : (byte)0;
|
||||
}
|
||||
|
||||
ulong pow = Pow2Mod(prime, divisor);
|
||||
hits[index] = pow == 1UL ? (byte)1 : (byte)0;
|
||||
}
|
||||
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
|
||||
{
|
||||
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
|
||||
{
|
||||
|
||||
Reference in new issue
Block a user