From d20362a2fc595a3445380b610d3c1f06e8c82ba8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Fri, 24 Oct 2025 09:20:17 +0200 Subject: [PATCH] Switch GPU by-divisor stepper to standard modular arithmetic --- PerfectNumbers.Core/DivisorCycleCache.cs | 4 +- .../Gpu/ExponentRemainderStepperGpu.cs | 190 ++++++++++++------ .../Gpu/Kernels/DivisorByDivisorKernels.cs | 164 +-------------- 3 files changed, 133 insertions(+), 225 deletions(-) diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index edc92d9..98fe8bf 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -79,7 +79,7 @@ public sealed class DivisorCycleCache throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } - ObserveCycleCacheDivisor(divisor); + // ObserveCycleCacheDivisor(divisor); return cached; } @@ -114,7 +114,7 @@ public sealed class DivisorCycleCache throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } - ObserveCycleCacheDivisor(divisor); + // ObserveCycleCacheDivisor(divisor); cycles[i] = cached; } else diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs index 9d652f5..0d74bc5 100644 --- a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs +++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs @@ -1,4 +1,6 @@ using System.Runtime.CompilerServices; +using ILGPU.Algorithms; +using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -8,22 +10,18 @@ namespace PerfectNumbers.Core.Gpu; /// internal struct ExponentRemainderStepperGpu { - private readonly MontgomeryDivisorData _divisor; + private const int WindowSizeMax = 8; + private readonly ulong _modulus; - private readonly ulong _nPrime; - private readonly ulong _montgomeryOne; private ulong _previousExponent; - private ulong _currentMontgomery; + private ulong _currentResidue; [MethodImpl(MethodImplOptions.AggressiveInlining)] public ExponentRemainderStepperGpu(in MontgomeryDivisorData divisor) { - _divisor = divisor; _modulus = divisor.Modulus; - _nPrime = divisor.NPrime; - _montgomeryOne = divisor.MontgomeryOne; _previousExponent = 0UL; - _currentMontgomery = divisor.MontgomeryOne; + _currentResidue = 1UL; } public ulong PreviousExponent @@ -35,18 +33,18 @@ internal struct ExponentRemainderStepperGpu [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool MatchesDivisor(in MontgomeryDivisorData divisor) { - return _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne; + return _modulus == divisor.Modulus; } /// - /// Clears the cached Montgomery residue so the stepper can be reused for another divisor. + /// Clears the cached residue so the stepper can be reused for another divisor. /// Callers must reinitialize the stepper before consuming any residues after a reset. /// [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Reset() { _previousExponent = 0UL; - _currentMontgomery = _montgomeryOne; + _currentResidue = 1UL; } /// @@ -58,20 +56,20 @@ internal struct ExponentRemainderStepperGpu [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong InitializeGpu(ulong exponent) { - InitializeGpuState(exponent); - return ReduceCurrent(); + InitializeStateGpu(exponent); + return _currentResidue; } /// /// Initializes the GPU stepping state and reports whether the seed exponent already yields unity. /// /// First exponent evaluated for the current divisor. Must match the divisor passed to the constructor. - /// when the seed exponent produces a Montgomery residue equal to one. + /// when the seed exponent produces a residue equal to one. [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool InitializeGpuIsUnity(ulong exponent) + public bool InitializeIsUnityGpu(ulong exponent) { - InitializeGpuState(exponent); - return _currentMontgomery == _montgomeryOne; + InitializeStateGpu(exponent); + return _currentResidue == 1UL; } /// @@ -83,74 +81,140 @@ internal struct ExponentRemainderStepperGpu [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong ComputeNextGpu(ulong exponent) { - AdvanceGpu(exponent); - return ReduceCurrent(); + AdvanceStateGpu(exponent); + return _currentResidue; } /// /// Advances the GPU stepping state to the next exponent and reports whether it evaluates to unity. - /// Callers must seed the stepper through or and then + /// Callers must seed the stepper through or and then /// pass exponents in strictly ascending order. /// /// Next exponent in the ascending sequence. /// when the canonical residue equals one. [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool ComputeNextGpuIsUnity(ulong exponent) + public bool ComputeNextIsUnityGpu(ulong exponent) { - AdvanceGpu(exponent); - return _currentMontgomery == _montgomeryOne; - } - - /// - /// Seeds the stepper with a precomputed Montgomery residue produced by the caller. - /// The provided exponent becomes the new baseline and all subsequent calls must continue with greater exponents. - /// - /// Exponent associated with . - /// Residue in Montgomery form produced externally for the configured divisor. - /// Outputs whether the supplied residue equals Montgomery one. - /// when the residue is accepted and cached. - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity) - { - _currentMontgomery = montgomeryResult; - _previousExponent = exponent; - isUnity = _currentMontgomery == _montgomeryOne; - return true; - } - - /// - /// Advances the stepper using a caller-supplied Montgomery delta that represents the gap to the next exponent. - /// Callers must seed the stepper first, reuse the same divisor metadata, and continue feeding exponents in ascending order so the delta remains valid. - /// - /// Next exponent in the ascending sequence. - /// Montgomery residue delta corresponding to the exponent gap. - /// Outputs whether the resulting residue equals Montgomery one. - /// when the advance succeeds. - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity) - { - _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime); - _previousExponent = exponent; - isUnity = _currentMontgomery == _montgomeryOne; - return true; + AdvanceStateGpu(exponent); + return _currentResidue == 1UL; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void InitializeGpuState(ulong exponent) + private void InitializeStateGpu(ulong exponent) { - _currentMontgomery = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, exponent); + _currentResidue = Pow2ModWindowed(exponent, _modulus); _previousExponent = exponent; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void AdvanceGpu(ulong exponent) + private void AdvanceStateGpu(ulong exponent) { ulong delta = exponent - _previousExponent; - ulong multiplier = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, delta); - _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); + ulong multiplier = Pow2ModWindowed(delta, _modulus); + _currentResidue = _currentResidue.MulMod64Gpu(multiplier, _modulus); _previousExponent = exponent; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime); + internal static ulong Pow2ModWindowed(ulong exponent, ulong modulus) + { + if (exponent == 0UL) + { + // The GPU by-divisor flow never routes modulus ≤ 1 here because divisors follow q = 2kp + 1 ≥ 3. + 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 = result.MulMod64Gpu(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 = result.MulMod64Gpu(result, modulus); + } + + if (!processWindow) + { + continue; + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + ulong multiplier = ComputeWindowedOddPower(windowValue, modulus); + result = result.MulMod64Gpu(multiplier, modulus); + index = windowStart - 1; + } + + return result; + } + + [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) + { + return baseValue; + } + + ulong result = baseValue; + ulong remaining = (windowValue - 1UL) >> 1; + ulong squareBase = baseValue.MulMod64Gpu(baseValue, modulus); + + while (remaining != 0UL) + { + result = ((remaining & 1UL) != 0UL) ? result.MulMod64Gpu(squareBase, modulus) : result; + remaining >>= 1; + squareBase = (remaining != 0UL) ? squareBase.MulMod64Gpu(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 2f835b5..88295f1 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -1,4 +1,3 @@ -using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; @@ -9,7 +8,6 @@ namespace PerfectNumbers.Core.Gpu; internal static class DivisorByDivisorKernels { - private const int WindowSizeMax = 8; public static void CheckKernel( Index1D index, ArrayView divisors, @@ -20,7 +18,7 @@ internal static class DivisorByDivisorKernels int globalIndex = index; ulong modulus = divisors[globalIndex]; ulong exponent = exponents[globalIndex]; - ulong result = Pow2ModWindowedKernel(exponent, modulus); + ulong result = ExponentRemainderStepperGpu.Pow2ModWindowed(exponent, modulus); byte hit = result == 1UL ? (byte)1 : (byte)0; hits[globalIndex] = hit; @@ -30,160 +28,6 @@ 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) - { - 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 ComputeWindowedOddPower(ulong windowValue, ulong modulus) - { - ulong baseValue = 2UL % modulus; - if (windowValue == 1UL) - { - return baseValue; - } - - ulong result = baseValue; - ulong remaining = (windowValue - 1UL) >> 1; - ulong squareBase = MultiplyMod(baseValue, baseValue, modulus); - - while (remaining != 0UL) - { - if ((remaining & 1UL) != 0UL) - { - result = MultiplyMod(result, squareBase, modulus); - } - - remaining >>= 1; - if (remaining == 0UL) - { - break; - } - - squareBase = MultiplyMod(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; - } - - [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); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong MultiplyMod(ulong left, ulong right, ulong modulus) - { - if (left == 0UL || right == 0UL) - { - return 0UL; - } - - ulong threshold = ulong.MaxValue / left; - if (right <= threshold) - { - ulong product = left * right; - return product % modulus; - } - - return MultiplyModWide(left, right, modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong MultiplyModWide(ulong left, ulong right, ulong modulus) - { - ulong result = 0UL; - ulong addend = left; - ulong remaining = right; - - while (remaining != 0UL) - { - if ((remaining & 1UL) != 0UL) - { - result += addend; - if (result >= modulus) - { - result -= modulus; - } - } - - addend <<= 1; - if (addend >= modulus) - { - addend -= modulus; - } - - remaining >>= 1; - } - - return result; - } - public static void EvaluateDivisorWithStepperKernel( Index1D _, MontgomeryDivisorData divisor, @@ -202,7 +46,7 @@ internal static class DivisorByDivisorKernels var exponentStepper = new ExponentRemainderStepperGpu(divisor); ulong firstPrime = exponents[0]; - bool firstUnity = exponentStepper.InitializeGpuIsUnity(firstPrime); + bool firstUnity = exponentStepper.InitializeIsUnityGpu(firstPrime); // EvenPerfectBitScanner only schedules divisors with non-zero cycle lengths (equal to the tested prime), // so the no-cycle path remains commented out to keep the kernel branch-free. // if (divisorCycle == 0UL) @@ -211,7 +55,7 @@ internal static class DivisorByDivisorKernels // for (int i = 1; i < length; i++) // { // ulong primeNoCycle = exponents[i]; - // bool unityNoCycle = exponentStepper.ComputeNextGpuIsUnity(primeNoCycle); + // bool unityNoCycle = exponentStepper.ComputeNextIsUnityGpu(primeNoCycle); // hits[i] = unityNoCycle ? (byte)1 : (byte)0; // } // return; @@ -234,7 +78,7 @@ internal static class DivisorByDivisorKernels cycleRemainder = remainderValue.Low; hits[i] = cycleRemainder == 0UL - ? (exponentStepper.ComputeNextGpuIsUnity(prime) ? (byte)1 : (byte)0) + ? (exponentStepper.ComputeNextIsUnityGpu(prime) ? (byte)1 : (byte)0) : (byte)0; } }