From 1916196e89711c85c5ce9e5b30e381a7d987b2aa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= <marek@lemanczyk.eu>
Date: Mon, 13 Oct 2025 14:21:28 +0200
Subject: [PATCH] Introduce windowed plain pow2mod helper for CPU checks

---
 PerfectNumbers.Core/MersenneDivisorCycles.cs |  4 +-
 PerfectNumbers.Core/ULongExtensions.cs       | 93 ++++++++++++++++++++
 2 files changed, 95 insertions(+), 2 deletions(-)

diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs
index b058e68..8266171 100644
--- a/PerfectNumbers.Core/MersenneDivisorCycles.cs
+++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs
@@ -109,7 +109,7 @@ public class MersenneDivisorCycles
             return cycleLength == exponent;
         }
 
-        return exponent.Pow2ModBinaryCpu(divisorData.Modulus) == 1UL;
+        return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL;
     }
 
     public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -133,7 +133,7 @@ public class MersenneDivisorCycles
             }
         }
 
-        return exponent.Pow2ModBinaryCpu(divisorData.Modulus) == 1UL;
+        return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL;
     }
 
     private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs
index 8fe7905..d8f074c 100644
--- a/PerfectNumbers.Core/ULongExtensions.cs
+++ b/PerfectNumbers.Core/ULongExtensions.cs
@@ -312,12 +312,81 @@ public static class ULongExtensions
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong Pow2ModBinaryCpu(this ulong exponent, ulong modulus)
+    {
+        return Pow2ModWindowedCpu(exponent, modulus);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus)
     {
         if (modulus <= 1UL)
         {
             return 0UL;
         }
 
+        if (exponent == 0UL)
+        {
+            return 1UL % modulus;
+        }
+
+        if (exponent <= Pow2WindowFallbackThreshold)
+        {
+            return Pow2ModBinaryFallback(exponent, modulus);
+        }
+
+        int bitLength = GetPortableBitLength(exponent);
+        int windowSize = GetWindowSize(bitLength);
+        int oddPowerCount = 1 << (windowSize - 1);
+
+        Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
+            ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]
+            : new ulong[oddPowerCount];
+        // Reuse oddPowerStorage to expose only the populated slice and avoid another local span.
+        oddPowerStorage = oddPowerStorage[..oddPowerCount];
+        InitializePlainOddPowers(modulus, oddPowerStorage);
+
+        ulong result = 1UL % modulus;
+        int index = bitLength - 1;
+
+        while (index >= 0)
+        {
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result = result.MulMod64(result, modulus);
+                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 = 0; square < windowLength; square++)
+            {
+                result = result.MulMod64(result, modulus);
+            }
+
+            ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL);
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            result = result.MulMod64(oddPowerStorage[tableIndex], modulus);
+
+            index = windowStart - 1;
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Pow2ModBinaryFallback(ulong exponent, ulong modulus)
+    {
         ulong remainingExponent = exponent;
         ulong result = 1UL % modulus;
         ulong baseValue = 2UL % modulus;
@@ -610,6 +679,30 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void InitializePlainOddPowers(ulong modulus, Span<ulong> destination)
+    {
+        int oddPowerCount = destination.Length;
+        if (oddPowerCount == 0)
+        {
+            return;
+        }
+
+        ulong baseValue = 2UL % modulus;
+        destination[0] = baseValue;
+        if (oddPowerCount == 1)
+        {
+            return;
+        }
+
+        ulong square = baseValue.MulMod64(baseValue, modulus);
+        for (int i = 1; i < oddPowerCount; i++)
+        {
+            // Reuse baseValue to accumulate successive odd powers and avoid allocating another local.
+            baseValue = baseValue.MulMod64(square, modulus);
+            destination[i] = baseValue;
+        }
+    }
+
     private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> destination)
     {
         int oddPowerCount = destination.Length;