diff --git a/EvenPerfectBitScanner.Benchmarks/GpuModularArithmeticBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuModularArithmeticBenchmarks.cs new file mode 100644 index 0000000..8b78308 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuModularArithmeticBenchmarks.cs @@ -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; + } +} diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 81b58e9..ed8ee1e 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -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 states = new(primes.Count); + List 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.Shared.Rent(stateCount); ulong[] primeBuffer = ArrayPool.Shared.Rent(stateCount); int[] indexBuffer = ArrayPool.Shared.Rent(stateCount); - List completions = new(stateCount); - List composites = new(stateCount); + PendingResult[] completionsBuffer = ArrayPool.Shared.Rent(stateCount); + PendingResult[] compositesBuffer = ArrayPool.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.Shared.Return(compositesBuffer, clearArray: true); + ArrayPool.Shared.Return(completionsBuffer, clearArray: true); ArrayPool.Shared.Return(indexBuffer, clearArray: true); ArrayPool.Shared.Return(primeBuffer, clearArray: true); ArrayPool.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(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 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 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; + } } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 604b0cf..a2da814 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -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, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelByPrimeCache = new(); + public readonly struct MontgomeryDivisorData + { + public readonly ulong Modulus; + public readonly ulong NPrime; + public readonly ulong MontgomeryOne; + public readonly ulong MontgomeryTwo; - private Action, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); - private Action, ArrayView> GetKernelByPrime(Accelerator accelerator) => - _kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernelByPrime)); + public MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo) + { + Modulus = modulus; + NPrime = nPrime; + MontgomeryOne = montgomeryOne; + MontgomeryTwo = montgomeryTwo; + } + } + + private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelByPrimeCache = new(); + + private Action, ArrayView> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); + private Action, ArrayView> GetKernelByPrime(Accelerator accelerator) => + _kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernelByPrime)); public int GpuBatchSize { @@ -90,15 +107,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester var kernel = GetKernel(accelerator); int batchCapacity = _gpuBatchSize; - var divisorsBuffer = accelerator.Allocate1D(batchCapacity); - var hitsBuffer = accelerator.Allocate1D(batchCapacity); + var divisorsBuffer = accelerator.Allocate1D(batchCapacity); + var hitsBuffer = accelerator.Allocate1D(batchCapacity); - ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); - byte[] hits = ArrayPool.Shared.Rent(batchCapacity); + ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); + byte[] hits = ArrayPool.Shared.Rent(batchCapacity); + MontgomeryDivisorData[] divisorData = ArrayPool.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 divisorDataSpan = divisorData.AsSpan(0, batchSize); + var divisorView = divisorsBuffer.View.SubView(0, batchSize); + var hitsView = hitsBuffer.View.SubView(0, batchSize); - Span divisorSpan = divisors.AsSpan(0, batchSize); - Span hitsSpan = hits.AsSpan(0, batchSize); + Span divisorSpan = divisors.AsSpan(0, batchSize); + Span 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.Shared.Return(hits, clearArray: true); - ArrayPool.Shared.Return(divisors, clearArray: true); - hitsBuffer.Dispose(); - divisorsBuffer.Dispose(); - gpuLease.Dispose(); - } + ArrayPool.Shared.Return(hits, clearArray: true); + ArrayPool.Shared.Return(divisors, clearArray: true); + ArrayPool.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 divisors, ArrayView hits) - { - ulong divisor = divisors[index]; - if (divisor <= 1UL || (divisor & 1UL) == 0UL) - { - hits[index] = 0; - return; - } + private static void CheckKernel(Index1D index, ulong prime, ArrayView divisors, ArrayView 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, ArrayView> _kernel; + private readonly Action, ArrayView> _kernel; private readonly MemoryBuffer1D _primesBuffer; private readonly MemoryBuffer1D _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(owner._gpuBatchSize); _hitsBuffer = _accelerator.Allocate1D(owner._gpuBatchSize); _primesHost = ArrayPool.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 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 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 primes, ArrayView 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 primes, ArrayView 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 {