diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 15d4023..db03894 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -1811,7 +1811,57 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { ulong exponent = exponents[index]; - remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false); + ulong modulus = divisor.Modulus; + ulong nPrime = divisor.NPrime; + ulong result = divisor.MontgomeryOne; + ulong baseValue = divisor.MontgomeryTwo; + ulong remainingExponent = exponent; + + if (remainingExponent == 0UL) + { + remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime); + return; + } + + while (remainingExponent != 0UL) + { + if ((remainingExponent & 1UL) != 0UL) + { + result = result.MontgomeryMultiply(baseValue, modulus, nPrime); + } + + remainingExponent >>= 1; + if (remainingExponent == 0UL) + { + break; + } + + // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term. + baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); + } + + remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime); + + // Plain modular exponentiation without Montgomery reduction kept for reference. + // ulong plainResult = 1UL % modulus; + // ulong plainBase = 2UL % modulus; + // ulong plainExponent = exponent; + // while (plainExponent != 0UL) + // { + // if ((plainExponent & 1UL) != 0UL) + // { + // plainResult = (plainResult * plainBase) % modulus; + // } + // + // plainExponent >>= 1; + // if (plainExponent == 0UL) + // { + // break; + // } + // + // plainBase = (plainBase * plainBase) % modulus; + // } + // remainders[index] = plainResult; } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator)