diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 2abfb65..7d5e482 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -1,3 +1,4 @@ +using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; @@ -8,21 +9,19 @@ namespace PerfectNumbers.Core.Gpu; internal static class DivisorByDivisorKernels { + private const int WindowSizeMax = 8; public static void CheckKernel( Index1D index, - ArrayView divisors, + ArrayView divisors, ArrayView exponents, ArrayView hits, ArrayView firstHit) { int globalIndex = index; - MontgomeryDivisorData divisor = divisors[globalIndex]; - ulong modulus = divisor.Modulus; - // The GPU by-divisor pipeline only materializes odd moduli greater than one (q = 2kp + 1), - // so the defensive guard for invalid values stays disabled here to keep the kernel branch-free. + ulong modulus = divisors[globalIndex]; ulong exponent = exponents[globalIndex]; - ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent); - byte hit = montgomeryResult == 1UL ? (byte)1 : (byte)0; + ulong result = Pow2ModWindowedKernel(exponent, modulus); + byte hit = result == 1UL ? (byte)1 : (byte)0; hits[globalIndex] = hit; if (hit != 0) @@ -31,6 +30,100 @@ internal static class DivisorByDivisorKernels } } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2ModWindowedKernel(ulong exponent, ulong modulus) + { + int bitLength = GetBitLength(exponent); + int windowSize = GetWindowSize(bitLength); + ulong result = 1UL; + int index = bitLength - 1; + + while (index >= 0) + { + result = ULongExtensions.MulMod64Gpu(result, result, modulus); + + if (((exponent >> index) & 1UL) == 0UL) + { + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowLength = index - windowStart + 1; + + for (int square = 1; square < windowLength; square++) + { + result = ULongExtensions.MulMod64Gpu(result, result, modulus); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + ulong multiplier = ComputeWindowedOddPower(windowValue, modulus); + result = ULongExtensions.MulMod64Gpu(result, multiplier, modulus); + + index = windowStart - 1; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeWindowedOddPower(ulong windowValue, ulong modulus) + { + ulong baseValue = 2UL % modulus; + ulong result = baseValue; + ulong remaining = (windowValue - 1UL) >> 1; + ulong squareBase = ULongExtensions.MulMod64Gpu(baseValue, baseValue, modulus); + + while (remaining != 0UL) + { + if ((remaining & 1UL) != 0UL) + { + result = ULongExtensions.MulMod64Gpu(result, squareBase, modulus); + } + + remaining >>= 1; + + if (remaining == 0UL) + { + break; + } + + squareBase = ULongExtensions.MulMod64Gpu(squareBase, squareBase, modulus); + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetBitLength(ulong value) + { + return 64 - XMath.LeadingZeroCount(value); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetWindowSize(int bitLength) + { + int window = WindowSizeMax; + window = bitLength <= 671 ? 7 : window; + window = bitLength <= 239 ? 6 : window; + window = bitLength <= 79 ? 5 : window; + window = bitLength <= 23 ? 4 : window; + int clamped = bitLength >= 1 ? bitLength : 1; + window = bitLength <= 6 ? clamped : window; + return window; + } public static void EvaluateDivisorWithStepperKernel( Index1D _, MontgomeryDivisorData divisor, diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 367be54..9c0f291 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -23,14 +23,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private ulong _divisorLimit; // private bool _isConfigured; - private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView>> _kernelCache = new(); + private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView>> _kernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> _stepperKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); - private Action, ArrayView, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(DivisorByDivisorKernels.CheckKernel)); + private Action, ArrayView, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(DivisorByDivisorKernels.CheckKernel)); private Action, ArrayView> GetStepperKernel(Accelerator accelerator) => _stepperKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.EvaluateDivisorWithStepperKernel)); @@ -112,7 +112,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.HitIndexBuffer, resources.Divisors, resources.Exponents, - resources.DivisorData, out lastProcessed, out coveredRange, out processedCount); @@ -163,14 +162,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static bool CheckDivisors( ulong prime, ulong allowedMax, - Action, ArrayView, ArrayView, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, + Action, ArrayView, ArrayView, ArrayView> kernel, + MemoryBuffer1D divisorsBuffer, MemoryBuffer1D exponentBuffer, MemoryBuffer1D hitsBuffer, MemoryBuffer1D hitIndexBuffer, in ulong[] divisors, ulong[] exponents, - MontgomeryDivisorData[] divisorData, out ulong lastProcessed, out bool coveredRange, out ulong processedCount) @@ -205,7 +203,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool ulongPool = ThreadStaticPools.UlongPool; ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity); Span filteredStorage = filteredDivisors.AsSpan(); - Span divisorDataStorage = divisorData.AsSpan(); Span divisorStorage = divisors.AsSpan(); Span exponentStorage = exponents.AsSpan(); @@ -254,7 +251,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Span filteredDivisorsSpan = filteredStorage[..filteredCount]; - Span divisorDataSpan = divisorDataStorage; Span divisorSpan = divisorStorage; Span exponentSpan = exponentStorage; @@ -271,7 +267,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv divisorSpan[admissibleCount] = divisorValue; exponentSpan[admissibleCount] = prime; - divisorDataSpan[admissibleCount] = divisorDataValue; admissibleCount++; } @@ -281,14 +276,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } Span activeDivisors = divisorSpan[..admissibleCount]; - Span activeDivisorData = divisorDataSpan[..admissibleCount]; Span activeExponents = exponentSpan[..admissibleCount]; - ArrayView1D divisorView = divisorsBuffer.View.SubView(0, admissibleCount); + ArrayView1D divisorView = divisorsBuffer.View.SubView(0, admissibleCount); ArrayView1D exponentView = exponentBuffer.View.SubView(0, admissibleCount); ArrayView1D hitsView = hitsBuffer.View.SubView(0, admissibleCount); ArrayView1D hitIndexView = hitIndexBuffer.View.SubView(0, 1); - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(activeDivisorData), admissibleCount); + divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(activeDivisors), admissibleCount); exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(activeExponents), admissibleCount); int sentinel = int.MaxValue; hitIndexView.SubView(0, 1).CopyFromCPU(ref sentinel, 1); @@ -586,18 +580,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { int actualCapacity = Math.Max(1, capacity); - DivisorsBuffer = accelerator.Allocate1D(actualCapacity); + DivisorsBuffer = accelerator.Allocate1D(actualCapacity); ExponentBuffer = accelerator.Allocate1D(actualCapacity); HitsBuffer = accelerator.Allocate1D(actualCapacity); HitIndexBuffer = accelerator.Allocate1D(1); ArrayPool ulongPool = ThreadStaticPools.UlongPool; Divisors = ulongPool.Rent(actualCapacity); Exponents = ulongPool.Rent(actualCapacity); - DivisorData = ThreadStaticPools.MontgomeryDivisorDataPool.Rent(actualCapacity); Capacity = actualCapacity; } - internal readonly MemoryBuffer1D DivisorsBuffer; + internal readonly MemoryBuffer1D DivisorsBuffer; internal readonly MemoryBuffer1D ExponentBuffer; @@ -609,8 +602,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly ulong[] Exponents; - internal readonly MontgomeryDivisorData[] DivisorData; - internal readonly int Capacity; public void Dispose() @@ -620,10 +611,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv HitsBuffer.Dispose(); HitIndexBuffer.Dispose(); ArrayPool ulongPool = ThreadStaticPools.UlongPool; - ArrayPool divisorPool = ThreadStaticPools.MontgomeryDivisorDataPool; ulongPool.Return(Divisors, clearArray: false); ulongPool.Return(Exponents, clearArray: false); - divisorPool.Return(DivisorData, clearArray: false); } }