diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs index c4c979f..dbfa2c1 100644 --- a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs +++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs @@ -1,5 +1,4 @@ using System.Runtime.CompilerServices; -using ILGPU.Algorithms; using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -10,18 +9,19 @@ namespace PerfectNumbers.Core.Gpu; /// internal struct ExponentRemainderStepperGpu { - private const int WindowSizeMax = 8; - private readonly ulong _modulus; + private readonly ReadOnlyGpuUInt128 _modulusWide; public ulong PreviousExponent; - private ulong _currentResidue; + private GpuUInt128 _currentResidue; [MethodImpl(MethodImplOptions.AggressiveInlining)] public ExponentRemainderStepperGpu(in MontgomeryDivisorData divisor) { - _modulus = divisor.Modulus; + ulong modulus = divisor.Modulus; + _modulus = modulus; + _modulusWide = new ReadOnlyGpuUInt128(modulus); PreviousExponent = 0UL; - _currentResidue = 1UL; + _currentResidue = new GpuUInt128(1UL); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -38,7 +38,7 @@ internal struct ExponentRemainderStepperGpu public void Reset() { PreviousExponent = 0UL; - _currentResidue = 1UL; + _currentResidue = new GpuUInt128(1UL); } /// @@ -51,7 +51,7 @@ internal struct ExponentRemainderStepperGpu public ulong InitializeGpu(ulong exponent) { InitializeStateGpu(exponent); - return _currentResidue; + return _currentResidue.Low; } /// @@ -63,7 +63,7 @@ internal struct ExponentRemainderStepperGpu public bool InitializeIsUnityGpu(ulong exponent) { InitializeStateGpu(exponent); - return _currentResidue == 1UL; + return _currentResidue.Low == 1UL; } /// @@ -76,7 +76,7 @@ internal struct ExponentRemainderStepperGpu public ulong ComputeNextGpu(ulong exponent) { AdvanceStateGpu(exponent); - return _currentResidue; + return _currentResidue.Low; } /// @@ -90,13 +90,14 @@ internal struct ExponentRemainderStepperGpu public bool ComputeNextIsUnityGpu(ulong exponent) { AdvanceStateGpu(exponent); - return _currentResidue == 1UL; + return _currentResidue.Low == 1UL; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private void InitializeStateGpu(ulong exponent) { - _currentResidue = Pow2ModWindowed(exponent, _modulus); + ulong residue = ULongExtensions.Pow2ModWindowedGpuKernel(exponent, _modulus); + _currentResidue = new GpuUInt128(residue); PreviousExponent = exponent; } @@ -104,120 +105,8 @@ internal struct ExponentRemainderStepperGpu private void AdvanceStateGpu(ulong exponent) { ulong delta = exponent - PreviousExponent; - ulong multiplier = Pow2ModWindowed(delta, _modulus); - _currentResidue = MultiplyMod(_currentResidue, multiplier, _modulus); + ulong multiplier = ULongExtensions.Pow2ModWindowedGpuKernel(delta, _modulus); + _currentResidue.MulMod(multiplier, _modulusWide); PreviousExponent = exponent; } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2ModWindowed(ulong exponent, ulong modulus) - { - // EvenPerfectBitScanner never provides zero exponents on the GPU path, so the guard stays commented out. - // if (exponent == 0UL) - // { - // return 1UL; - // } - - int bitLength = GetBitLength(exponent); - int windowSize = GetWindowSize(bitLength); - ulong result = 1UL; - int index = bitLength - 1; - - while (index >= 0) - { - ulong currentBit = (exponent >> index) & 1UL; - ulong squared = MultiplyMod(result, result, modulus); - bool processWindow = currentBit != 0UL; - - result = processWindow ? result : squared; - index -= (int)(currentBit ^ 1UL); - - int windowStartCandidate = index - windowSize + 1; - int negativeMask = windowStartCandidate >> 31; - windowStartCandidate &= ~negativeMask; - - int windowStart = processWindow ? GetNextSetBitIndex(exponent, windowStartCandidate) : windowStartCandidate; - int windowLength = processWindow ? index - windowStart + 1 : 0; - - for (int square = 0; square < windowLength; square++) - { - result = MultiplyMod(result, result, modulus); - } - - if (!processWindow) - { - continue; - } - - ulong mask = (1UL << windowLength) - 1UL; - ulong windowValue = (exponent >> windowStart) & mask; - ulong multiplier = ComputeWindowedOddPower(windowValue, modulus); - result = MultiplyMod(result, multiplier, modulus); - index = windowStart - 1; - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong MultiplyMod(ulong left, ulong right, ulong modulus) - { - // GpuUInt128.MulMod is the fastest GPU-compatible multiply-reduce helper available to this kernel path. - var product = new GpuUInt128(left); - return product.MulMod(right, modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeWindowedOddPower(ulong windowValue, ulong modulus) - { - // Production divisors on the GPU path always satisfy modulus ≥ 3 (q = 2kp + 1), so the base stays 2 without a reduction. - ulong baseValue = 2UL; - if (windowValue == 1UL) - { - // Single-bit windows occur frequently in production primes, so keep the fast return for window value 1. - return baseValue; - } - - ulong result = baseValue; - ulong remaining = (windowValue - 1UL) >> 1; - ulong squareBase = MultiplyMod(baseValue, baseValue, modulus); - - while (remaining != 0UL) - { - result = ((remaining & 1UL) != 0UL) ? MultiplyMod(result, squareBase, modulus) : result; - remaining >>= 1; - squareBase = (remaining != 0UL) ? MultiplyMod(squareBase, squareBase, modulus) : squareBase; - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetBitLength(ulong value) - { - return value == 0UL ? 0 : 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; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetNextSetBitIndex(ulong exponent, int startIndex) - { - ulong guard = (ulong)(((long)startIndex - 64) >> 63); - int shift = startIndex & 63; - ulong mask = (~0UL << shift) & guard; - ulong masked = exponent & mask; - return XMath.TrailingZeroCount(masked); - } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 88295f1..4fd04c9 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -18,7 +18,7 @@ internal static class DivisorByDivisorKernels int globalIndex = index; ulong modulus = divisors[globalIndex]; ulong exponent = exponents[globalIndex]; - ulong result = ExponentRemainderStepperGpu.Pow2ModWindowed(exponent, modulus); + ulong result = ULongExtensions.Pow2ModWindowedGpuKernel(exponent, modulus); byte hit = result == 1UL ? (byte)1 : (byte)0; hits[globalIndex] = hit; diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index d831097..c1d193c 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -108,6 +108,77 @@ public static partial class ULongExtensions return result; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong Pow2ModWindowedGpuKernel(ulong exponent, ulong modulus) + { + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + ReadOnlyGpuUInt128 modulusWide = new(modulus); + GpuUInt128 result = new(1UL); + int index = bitLength - 1; + + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + if (currentBit == 0UL) + { + result.MulMod(result.AsReadOnly(), modulusWide); + index--; + continue; + } + + int windowStartCandidate = index - windowSize + 1; + int negativeMask = windowStartCandidate >> 31; + windowStartCandidate &= ~negativeMask; + int windowStart = GetNextSetBitIndexGpu(exponent, windowStartCandidate); + int windowLength = index - windowStart + 1; + + for (int square = 0; square < windowLength; square++) + { + result.MulMod(result.AsReadOnly(), modulusWide); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + ulong multiplier = ComputeWindowedOddPowerGpuKernel(windowValue, modulusWide); + result.MulMod(multiplier, modulusWide); + + index = windowStart - 1; + } + + return result.Low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeWindowedOddPowerGpuKernel(ulong windowValue, in ReadOnlyGpuUInt128 modulus) + { + GpuUInt128 result = new(2UL); + if (windowValue == 1UL) + { + return 2UL; + } + + ulong remaining = (windowValue - 1UL) >> 1; + GpuUInt128 squareBase = new(2UL); + squareBase.MulMod(squareBase.AsReadOnly(), modulus); + + while (remaining != 0UL) + { + if ((remaining & 1UL) != 0UL) + { + result.MulMod(squareBase.AsReadOnly(), modulus); + } + + remaining >>= 1; + if (remaining != 0UL) + { + squareBase.MulMod(squareBase.AsReadOnly(), modulus); + } + } + + return result.Low; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorData divisor, ulong exponent) {