diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 3986136..44c9ebc 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -4,6 +4,7 @@ using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -24,14 +25,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester public readonly ulong MontgomeryTwo = montgomeryTwo; } - private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelByPrimeCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelByPrimeCache = new(); + private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _sessionPool = 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)); + 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 { @@ -86,37 +88,32 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); var accelerator = gpuLease.Accelerator; - var kernel = GetKernel(accelerator); + var kernel = GetKernel(accelerator); + BatchResources resources = RentBatchResources(accelerator, batchCapacity); - // TODO: Can we move the creation of the buffers to one of the parent methods, outside of all loops and initialize them once for all divisors / candidates? Instead of creating buffers for each thread separately, let's pre-allocate them and re-use in threads. - var divisorsBuffer = accelerator.Allocate1D(batchCapacity); - var hitsBuffer = accelerator.Allocate1D(batchCapacity); - ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); - byte[] hits = ArrayPool.Shared.Rent(batchCapacity); - MontgomeryDivisorData[] divisorData = ArrayPool.Shared.Rent(batchCapacity); - - composite = CheckDivisors( - prime, - allowedMax, - useCycles, - accelerator, - kernel, - divisorsBuffer, - hitsBuffer, - divisors, - hits, - divisorData, - out lastProcessed, - out coveredRange, - out processedCount - ); - - ArrayPool.Shared.Return(hits, clearArray: true); - ArrayPool.Shared.Return(divisors, clearArray: true); - ArrayPool.Shared.Return(divisorData, clearArray: true); - hitsBuffer.Dispose(); - divisorsBuffer.Dispose(); - gpuLease.Dispose(); + try + { + composite = CheckDivisors( + prime, + allowedMax, + useCycles, + accelerator, + kernel, + resources.DivisorsBuffer, + resources.HitsBuffer, + resources.Divisors, + resources.Hits, + resources.DivisorData, + out lastProcessed, + out coveredRange, + out processedCount + ); + } + finally + { + ReturnBatchResources(accelerator, resources); + gpuLease.Dispose(); + } if (processedCount > 0UL) { @@ -141,7 +138,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester ulong allowedMax, bool useCycles, Accelerator accelerator, - Action, ArrayView> kernel, + Action, ArrayView> kernel, MemoryBuffer1D divisorsBuffer, MemoryBuffer1D hitsBuffer, ulong[] divisors, @@ -164,18 +161,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester bool composite = false; bool processedAll = false; - ulong cycle, currentDivisor, divisor = 3UL, next, nextDivisor; + ulong currentDivisor, divisor = 3UL, next, nextDivisor; int batchSize, i; bool reachedEndInBatch; Span divisorDataSpan; Span divisorSpan; Span hitsSpan; - MersenneDivisorCycles divisorCycles = MersenneDivisorCycles.Shared; - ArrayView1D divisorView; - ArrayView1D hitsView; + ArrayView1D divisorView; + ArrayView1D hitsView; - while (divisor <= allowedMax) + while (divisor <= allowedMax) { batchSize = 0; reachedEndInBatch = false; @@ -184,34 +180,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester { currentDivisor = divisor; nextDivisor = currentDivisor + 2UL; - processedCount++; - - // TODO: Let's move all the pre-checks to GPU. - if (useCycles) - { - // TODO: Let's not use the divisor cache but instead calculate it on GPU. Let's move the entire check of cycle's divisibility to GPU. We should already be able to calculate in on GPU. There should be classes for that. - cycle = divisorCycles!.GetCycle(currentDivisor); - if (cycle == 0UL || prime % cycle != 0UL) - { - lastProcessed = currentDivisor; - - if (nextDivisor <= currentDivisor) - { - reachedEndInBatch = true; - break; - } - - if (nextDivisor > allowedMax) - { - divisor = nextDivisor; - reachedEndInBatch = true; - break; - } - - divisor = nextDivisor; - continue; - } - } + processedCount++; divisors[batchSize++] = currentDivisor; lastProcessed = currentDivisor; @@ -257,7 +226,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester } divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); - kernel(batchSize, prime, divisorView, hitsView); + kernel(batchSize, prime, useCycles ? (byte)1 : (byte)0, divisorView, hitsView); accelerator.Synchronize(); hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); @@ -336,42 +305,44 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester return; } - ulong total = _lastStatusDivisor + processedCount; - ulong messages = total / interval; - _lastStatusDivisor = total % interval; + ulong total = _lastStatusDivisor + processedCount; + _lastStatusDivisor = total % interval; - // TODO: Why do we output so many the same lines? Let's remove this printing to console. - for (ulong i = 0; i < messages; i++) - { - Console.WriteLine($"...processed by-divisor candidate = {lastProcessed}"); - } - } + // Removed noisy console output; status is tracked internally only. + } - private static void CheckKernel(Index1D index, ulong prime, ArrayView divisors, ArrayView hits) - { - MontgomeryDivisorData divisor = divisors[index]; - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } + private static void CheckKernel(Index1D index, ulong prime, byte useCycleChecks, ArrayView divisors, ArrayView hits) + { + MontgomeryDivisorData divisor = divisors[index]; + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits[index] = 0; + return; + } - hits[index] = Pow2Mod(prime, divisor) == 1UL ? (byte)1 : (byte)0; - } + if (useCycleChecks != 0) + { + ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); + if (cycle == 0UL || prime % cycle != 0UL) + { + hits[index] = 0; + return; + } + } + + hits[index] = prime.Pow2MontgomeryMod(divisor) == 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; - private readonly byte[] _hitsHost; - private readonly int[] _indexHost; - private bool _disposed; + private readonly MemoryBuffer1D _hitsBuffer; + private bool _disposed; internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) { @@ -381,10 +352,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester _kernel = owner.GetKernelByPrime(_accelerator); _primesBuffer = _accelerator.Allocate1D(owner._gpuBatchSize); _hitsBuffer = _accelerator.Allocate1D(owner._gpuBatchSize); - _primesHost = ArrayPool.Shared.Rent(owner._gpuBatchSize); - _hitsHost = ArrayPool.Shared.Rent(owner._gpuBatchSize); - _indexHost = ArrayPool.Shared.Rent(owner._gpuBatchSize); - } + } internal void Reset() { @@ -402,90 +370,38 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester MersenneNumberDivisorByDivisorGpuTester owner = _owner; int gpuBatchSize = owner._gpuBatchSize; bool useCycles = owner._useDivisorCycles; - byte[] hitsHost = _hitsHost; - int[] indexHost = _indexHost; - ulong[] primesHost = _primesHost; - Accelerator accelerator = _accelerator; - MemoryBuffer1D hitsBuffer = _hitsBuffer; - MemoryBuffer1D primesBuffer = _primesBuffer; - ArrayView1D primesBufferView = primesBuffer.View; - ArrayView1D hitsBufferView = hitsBuffer.View; - Action, ArrayView> kernel = _kernel; - ArrayView1D primesView; - ArrayView1D hitsView; - MersenneDivisorCycles? divisorCycles = MersenneDivisorCycles.Shared; - ulong cycle = useCycles ? divisorCycles!.GetCycle(divisor) : 0UL, primeValue; - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); - Span primeSpan = []; - ReadOnlySpan primesSlice; - Span indexSpan = []; - Span hitSpan = []; - Span hitsSlice; - int batchSize, currentSpanSize = -1, gpuCount, i, primesLength = primes.Length; + Accelerator accelerator = _accelerator; + MemoryBuffer1D hitsBuffer = _hitsBuffer; + MemoryBuffer1D primesBuffer = _primesBuffer; + ArrayView1D primesBufferView = primesBuffer.View; + ArrayView1D hitsBufferView = hitsBuffer.View; + Action, ArrayView> kernel = _kernel; + ArrayView1D primesView; + ArrayView1D hitsView; + byte useCyclesFlag = useCycles ? (byte)1 : (byte)0; + MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + ReadOnlySpan primesSlice; + Span hitsSlice; + int batchSize, primesLength = primes.Length; - while (offset < primesLength) - { - batchSize = Math.Min(gpuBatchSize, primesLength - offset); - primesSlice = primes.Slice(offset, batchSize); + while (offset < primesLength) + { + batchSize = Math.Min(gpuBatchSize, primesLength - offset); + primesSlice = primes.Slice(offset, batchSize); hitsSlice = hits.Slice(offset, batchSize); - if (useCycles) - { - if (batchSize != currentSpanSize) - { - primeSpan = primesHost.AsSpan(0, batchSize); - indexSpan = indexHost.AsSpan(0, batchSize); - hitSpan = hitsHost.AsSpan(0, batchSize); - currentSpanSize = batchSize; - } + primesView = primesBufferView.SubView(0, batchSize); + hitsView = hitsBufferView.SubView(0, batchSize); - gpuCount = 0; + primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); + kernel(batchSize, divisorData, useCyclesFlag, primesView, hitsView); + accelerator.Synchronize(); - for (i = 0; i < batchSize; i++) - { - primeValue = primesSlice[i]; - if (cycle == 0UL || primeValue % cycle != 0UL) - { - hitsSlice[i] = 0; - continue; - } + hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSlice), batchSize); - indexSpan[gpuCount] = i; - primeSpan[gpuCount] = primeValue; - gpuCount++; - } - - if (gpuCount > 0) - { - primesView = primesBufferView.SubView(0, gpuCount); - hitsView = hitsBufferView.SubView(0, gpuCount); - - primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primeSpan), gpuCount); - kernel(gpuCount, divisorData, primesView, hitsView); - accelerator.Synchronize(); - - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitSpan), gpuCount); - for (i = 0; i < gpuCount; i++) - { - hitsSlice[indexSpan[i]] = hitSpan[i]; - } - } - } - else - { - primesView = primesBufferView.SubView(0, batchSize); - hitsView = hitsBufferView.SubView(0, batchSize); - - primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); - kernel(batchSize, divisorData, primesView, hitsView); - accelerator.Synchronize(); - - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSlice), batchSize); - } - - offset += batchSize; - } - } + offset += batchSize; + } + } public void Dispose() { @@ -504,67 +420,29 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester _sessionPool.Add(session); } - // TODO: This method with all related should be moved to ULongExtensions and converted to extension method. We should already have similar methods. How this one is different? - private static ulong Pow2Mod(ulong exponent, in MontgomeryDivisorData divisor) - { - // TODO: Do we use the dependency and the residue for the previous p in the calculations? Let's make use of the previous residue for the previous p when calculating the residue for the same divisor but different candidate. The candidates should be processed in increasing order. - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return 0UL; - } + private static void CheckKernelByPrime(Index1D index, MontgomeryDivisorData divisor, byte useCycleChecks, ArrayView primes, ArrayView hits) + { + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits[index] = 0; + return; + } - ulong result = divisor.MontgomeryOne; - ulong baseVal = divisor.MontgomeryTwo; - ulong nPrime = divisor.NPrime; - while (exponent > 0UL) - { - if ((exponent & 1UL) != 0UL) - { - result = MontgomeryMultiply(result, baseVal, modulus, nPrime); - } + if (useCycleChecks != 0) + { + ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); + if (cycle == 0UL || primes[index] % cycle != 0UL) + { + hits[index] = 0; + return; + } + } - exponent >>= 1; - if (exponent == 0UL) - { - break; - } + hits[index] = primes[index].Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + } - baseVal = MontgomeryMultiply(baseVal, baseVal, modulus, nPrime); - } - - return MontgomeryMultiply(result, 1UL, modulus, nPrime); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong MontgomeryMultiply(ulong a, ulong b, ulong modulus, ulong nPrime) - { - ulong tLow = unchecked(a * b); - ulong m = unchecked(tLow * nPrime); - - // We're reusing tLow variable as the result to perform just a little bit better - tLow = unchecked(MultiplyHigh(a, b) + MultiplyHigh(m, modulus) + (unchecked(tLow + m * modulus) < tLow ? 1UL : 0UL)); - if (tLow >= modulus) - { - tLow -= modulus; - } - - return tLow; - } - - private static void CheckKernelByPrime(Index1D index, MontgomeryDivisorData divisor, ArrayView primes, ArrayView hits) - { - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - - hits[index] = Pow2Mod(primes[index], divisor) == 1UL ? (byte)1 : (byte)0; - } - - private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) + private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) { if (modulus <= 1UL || (modulus & 1UL) == 0UL) { @@ -593,35 +471,112 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester 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 mid = xHigh * yLow + ((xLow * yLow) >> 32); - return xHigh * yHigh + (mid >> 32) + (((mid & 0xFFFFFFFFUL) + xLow * yHigh) >> 32); - } - - public ulong DivisorLimit - { - get + public ulong DivisorLimit + { + get { lock (_sync) { if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } - return _divisorLimit; - } - } - } + return _divisorLimit; + } + } + } - public ulong GetAllowedMaxDivisor(ulong prime) + private BatchResources RentBatchResources(Accelerator accelerator, int capacity) + { + var bag = _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()); + if (bag.TryTake(out BatchResources? resources)) + { + resources.EnsureCapacity(capacity); + return resources; + } + + return new BatchResources(accelerator, capacity); + } + + private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) + { + _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()).Add(resources); + } + + private sealed class BatchResources : IDisposable + { + private readonly Accelerator _accelerator; + + internal BatchResources(Accelerator accelerator, int capacity) + { + _accelerator = accelerator; + DivisorsBuffer = accelerator.Allocate1D(capacity); + HitsBuffer = accelerator.Allocate1D(capacity); + Divisors = ArrayPool.Shared.Rent(capacity); + Hits = ArrayPool.Shared.Rent(capacity); + DivisorData = ArrayPool.Shared.Rent(capacity); + Capacity = capacity; + } + + internal MemoryBuffer1D DivisorsBuffer { get; private set; } + + internal MemoryBuffer1D HitsBuffer { get; private set; } + + internal ulong[] Divisors { get; private set; } + + internal byte[] Hits { get; private set; } + + internal MontgomeryDivisorData[] DivisorData { get; private set; } + + internal int Capacity { get; private set; } + + public void EnsureCapacity(int requiredCapacity) + { + if (requiredCapacity <= Capacity) + { + return; + } + + Resize(requiredCapacity); + } + + private void Resize(int newCapacity) + { + DivisorsBuffer.Dispose(); + HitsBuffer.Dispose(); + ArrayPool.Shared.Return(Divisors, clearArray: false); + ArrayPool.Shared.Return(Hits, clearArray: false); + ArrayPool.Shared.Return(DivisorData, clearArray: false); + + DivisorsBuffer = _accelerator.Allocate1D(newCapacity); + HitsBuffer = _accelerator.Allocate1D(newCapacity); + Divisors = ArrayPool.Shared.Rent(newCapacity); + Hits = ArrayPool.Shared.Rent(newCapacity); + DivisorData = ArrayPool.Shared.Rent(newCapacity); + Capacity = newCapacity; + } + + public void Dispose() + { + DivisorsBuffer.Dispose(); + HitsBuffer.Dispose(); + ArrayPool.Shared.Return(Divisors, clearArray: false); + ArrayPool.Shared.Return(Hits, clearArray: false); + ArrayPool.Shared.Return(DivisorData, clearArray: false); + } + } + + private sealed class AcceleratorReferenceComparer : IEqualityComparer + { + public static AcceleratorReferenceComparer Instance { get; } = new(); + + public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y); + + public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); + } + + public ulong GetAllowedMaxDivisor(ulong prime) { lock (_sync) { diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index cb89e14..726dc35 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -1,5 +1,6 @@ using System.Numerics; using System.Runtime.CompilerServices; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -157,11 +158,11 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b); - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulHigh(this ulong x, ulong y) - { - ulong xLow = (uint)x; - ulong xHigh = x >> 32; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulHigh(this ulong x, ulong y) + { + ulong xLow = (uint)x; + ulong xHigh = x >> 32; ulong yLow = (uint)y; ulong yHigh = y >> 32; @@ -174,11 +175,11 @@ public static class ULongExtensions (uint)w1 + (uint)w2 ) >> 32; - } + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) - { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) + { ulong result = 1UL; value %= modulus; @@ -196,8 +197,57 @@ public static class ULongExtensions return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) + { + ulong tLow = unchecked(a * b); + ulong m = unchecked(tLow * nPrime); + ulong mTimesModulusLow = unchecked(m * modulus); + + ulong result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL)); + if (result >= modulus) + { + result -= modulus; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryMod(this ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + { + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return 0UL; + } + + ulong result = divisor.MontgomeryOne; + ulong baseVal = divisor.MontgomeryTwo; + ulong nPrime = divisor.NPrime; + ulong remainingExponent = exponent; + + while (remainingExponent > 0UL) + { + if ((remainingExponent & 1UL) != 0UL) + { + result = result.MontgomeryMultiply(baseVal, modulus, nPrime); + } + + remainingExponent >>= 1; + if (remainingExponent == 0UL) + { + break; + } + + baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); + } + + return result.MontgomeryMultiply(1UL, modulus, nPrime); + } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus)