Merge pull request #175 from mlemanczyk/codex/optimize-evenperfectbitscanner-for-reduced-gc

Reduce GC pressure in divisor checking pipeline

EvenPerfectBitScanner/Program.cs

@@ -1,5 +1,6 @@
 using System.Buffers;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Globalization;
 using System.Numerics;
 using System.Runtime.CompilerServices;
@@ -45,6 +46,11 @@ internal static class Program
     private static readonly Optimized PrimeIterator = new();
     private static ulong _byDivisorStartPrime;
 
+    private static bool _testMode;
+    private static int _testTargetPrimeCount;
+    private static int _testProcessedPrimeCount;
+    private const string TestTimeFileName = "even_perfect_test_time.txt";
+
     [ThreadStatic]
     private static bool _lastCompositeP;
 
@@ -95,6 +101,10 @@ internal static class Program
         ulong remainder = 0UL;
         bool startPrimeProvided = false;
 
+        _testMode = false;
+        _testTargetPrimeCount = 0;
+        _testProcessedPrimeCount = 0;
+
         // NTT backend selection (GPU): reference vs staged
         for (; argIndex < args.Length; argIndex++)
         {
@@ -201,6 +211,10 @@ internal static class Program
             {
                 useModuloWorkaround = true;
             }
+            else if (arg.Equals("--test", StringComparison.OrdinalIgnoreCase))
+            {
+                _testMode = true;
+            }
             else if (arg.Equals("--use-order", StringComparison.OrdinalIgnoreCase))
             {
                 useOrder = true;
@@ -393,6 +407,19 @@ internal static class Program
             }
         }
 
+        int testPrimeCandidateLimit = 0;
+
+        if (_testMode)
+        {
+            currentP = 31UL;
+            remainder = currentP % 6UL;
+            filterFile = string.Empty;
+            maxPrimeLimit = UInt128.MaxValue;
+            maxPrimeConfigured = false;
+            startPrimeProvided = true;
+            testPrimeCandidateLimit = Math.Max(1, threadCount) * 3;
+        }
+
         if (useByDivisor)
         {
             _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL;
@@ -402,7 +429,7 @@ internal static class Program
             _byDivisorStartPrime = 0UL;
         }
 
-        if (useByDivisor && string.IsNullOrEmpty(filterFile))
+        if (useByDivisor && !_testMode && string.IsNullOrEmpty(filterFile))
         {
             Console.WriteLine("--mersenne=bydivisor requires --filter-p=<path>.");
             return;
@@ -620,17 +647,22 @@ internal static class Program
             File.WriteAllText(ResultsFileName, $"p,searchedMersenne,detailedCheck,passedAllTests{Environment.NewLine}");
         }
 
-        bool useFilter = !string.IsNullOrEmpty(filterFile);
+        bool useFilter = !_testMode && !string.IsNullOrEmpty(filterFile);
         HashSet<ulong> filter = [];
         List<ulong> byDivisorCandidates = [];
         ulong maxP = 0UL;
         ulong[] localFilter = Array.Empty<ulong>();
         int filterCount = 0;
-        if (useFilter)
+
+        if (_useByDivisorMode)
         {
-            Console.WriteLine("Loading filter...");
+            if (_testMode)
-            if (_useByDivisorMode)
             {
+                byDivisorCandidates = BuildTestPrimeCandidates(testPrimeCandidateLimit);
+            }
+            else
+            {
+                Console.WriteLine("Loading filter...");
                 int skippedByLimit = 0;
                 byDivisorCandidates = LoadByDivisorCandidates(filterFile, maxPrimeLimit, maxPrimeConfigured, out skippedByLimit);
                 if (maxPrimeConfigured)
@@ -647,63 +679,70 @@ internal static class Program
                     }
                 }
             }
-            else
+        }
+        else if (useFilter)
+        {
+            Console.WriteLine("Loading filter...");
+            // TODO: Rent this filter buffer from ArrayPool<ulong> and reuse it across reload batches so
+            // we do not allocate fresh arrays while replaying large result filters.
+            localFilter = new ulong[1024];
+            filterCount = 0;
+            bool addedWithinLimit = false;
+            int skippedByLimit = 0;
+            LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) =>
             {
-                // TODO: Rent this filter buffer from ArrayPool<ulong> and reuse it across reload batches so
+                if (!passedAllTests)
-                // we do not allocate fresh arrays while replaying large result filters.
-                localFilter = new ulong[1024];
-                filterCount = 0;
-                bool addedWithinLimit = false;
-                int skippedByLimit = 0;
-                LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) =>
                 {
-                    if (!passedAllTests)
+                    return;
-                    {
+                }
-                        return;
-                    }
 
-                    if (maxPrimeConfigured && (UInt128)p > maxPrimeLimit)
+                if (maxPrimeConfigured && (UInt128)p > maxPrimeLimit)
-                    {
+                {
-                        skippedByLimit++;
+                    skippedByLimit++;
-                        return;
+                    return;
-                    }
+                }
 
-                    localFilter[filterCount++] = p;
+                localFilter[filterCount++] = p;
-                    addedWithinLimit = true;
+                addedWithinLimit = true;
-                    if (p > maxP)
+                if (p > maxP)
-                    {
+                {
-                        maxP = p;
+                    maxP = p;
-                    }
+                }
 
-                    if (filterCount == 1024)
+                if (filterCount == 1024)
-                    {
-                        filter.AddRange(localFilter[..filterCount]);
-                        filterCount = 0;
-                        Console.WriteLine($"Added {p}");
-                    }
-                });
-
-                if (filterCount > 0)
                 {
                     filter.AddRange(localFilter[..filterCount]);
+                    filterCount = 0;
+                    Console.WriteLine($"Added {p}");
+                }
+            });
+
+            if (filterCount > 0)
+            {
+                filter.AddRange(localFilter[..filterCount]);
+            }
+
+            if (maxPrimeConfigured)
+            {
+                if (!addedWithinLimit)
+                {
+                    Console.WriteLine("No filter candidates fall within the --max-prime limit.");
+                    return;
                 }
 
-                if (maxPrimeConfigured)
+                if (skippedByLimit > 0)
                 {
-                    if (!addedWithinLimit)
+                    Console.WriteLine($"Skipped {skippedByLimit} filter candidate(s) above --max-prime limit.");
-                    {
-                        Console.WriteLine("No filter candidates fall within the --max-prime limit.");
-                        return;
-                    }
-
-                    if (skippedByLimit > 0)
-                    {
-                        Console.WriteLine($"Skipped {skippedByLimit} filter candidate(s) above --max-prime limit.");
-                    }
                 }
             }
         }
 
+
+        if (_testMode)
+        {
+            _testTargetPrimeCount = testPrimeCandidateLimit;
+        }
+
         if (_primeCount >= 2)
         {
             _primeFoundAfterInit = true;
@@ -716,6 +755,12 @@ internal static class Program
         // Configure batch size for GPU primality sieve
         PrimeTester.GpuBatchSize = gpuPrimeBatch;
 
+        Stopwatch? stopwatch = null;
+        if (_testMode)
+        {
+            stopwatch = Stopwatch.StartNew();
+        }
+
         Console.WriteLine("Starting scan...");
 
         if (_useByDivisorMode)
@@ -723,6 +768,12 @@ internal static class Program
             RunByDivisorMode(byDivisorCandidates, threadCount);
             FlushBuffer();
             StringBuilderPool.Return(_outputBuilder!);
+            if (stopwatch is not null)
+            {
+                stopwatch.Stop();
+                ReportTestTime(stopwatch.Elapsed);
+            }
+
             return;
         }
 
@@ -737,7 +788,14 @@ internal static class Program
             }
         }
 
+        if (_testMode)
+        {
+            testPrimeCandidateLimit = Math.Max(1, threadCount) * 3;
+            _testTargetPrimeCount = testPrimeCandidateLimit;
+        }
+
         _state = ((long)currentP << 3) | (long)remainder;
+
         Task[] tasks = new Task[threadCount];
 
         int taskIndex = 0;
@@ -807,8 +865,65 @@ internal static class Program
         }
 
         Task.WaitAll(tasks);
+
+        TimeSpan testElapsed = default;
+        bool reportTestTime = false;
+        if (stopwatch is not null)
+        {
+            stopwatch.Stop();
+            testElapsed = stopwatch.Elapsed;
+            reportTestTime = true;
+        }
+
         FlushBuffer();
         StringBuilderPool.Return(_outputBuilder!);
+
+        if (reportTestTime)
+        {
+            ReportTestTime(testElapsed);
+        }
+    }
+
+    private static List<ulong> BuildTestPrimeCandidates(int targetCount)
+    {
+        if (targetCount <= 0)
+        {
+            return [];
+        }
+
+        List<ulong> candidates = new(targetCount);
+        ulong previous = 29UL;
+        int remaining = targetCount;
+        ulong nextPrime = 0UL;
+
+        while (remaining > 0)
+        {
+            try
+            {
+                nextPrime = PrimeIterator.Next(in previous);
+            }
+            catch (InvalidOperationException)
+            {
+                break;
+            }
+
+            if (nextPrime < 31UL)
+            {
+                previous = nextPrime;
+                continue;
+            }
+
+            candidates.Add(nextPrime);
+            remaining--;
+            previous = nextPrime;
+        }
+
+        if (remaining > 0)
+        {
+            Console.WriteLine("Unable to populate the requested number of test primes before reaching the 64-bit limit.");
+        }
+
+        return candidates;
     }
 
     private static List<ulong> LoadByDivisorCandidates(string candidateFile, UInt128 maxPrimeLimit, bool maxPrimeConfigured, out int skippedByLimit)
@@ -1053,6 +1168,21 @@ internal static class Program
         }
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void OnPrimeCandidateConfirmed()
+    {
+        if (!_testMode)
+        {
+            return;
+        }
+
+        int processed = Interlocked.Increment(ref _testProcessedPrimeCount);
+        if (processed == _testTargetPrimeCount)
+        {
+            Volatile.Write(ref _limitReached, true);
+        }
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static bool TryAcquireConsoleSlot(bool primeCandidate, out bool primeFlag)
     {
@@ -1189,6 +1319,75 @@ internal static class Program
         }
     }
 
+    private static void ReportTestTime(TimeSpan elapsed)
+    {
+        double elapsedSeconds = elapsed.TotalSeconds;
+        double elapsedMilliseconds = elapsed.TotalMilliseconds;
+        Console.WriteLine($"Test elapsed time: {elapsedSeconds:F3} s");
+
+        string filePath = GetTestTimeFilePath();
+        bool hasPrevious = TryReadTestTime(filePath, out double previousMilliseconds);
+        string elapsedText = elapsedMilliseconds.ToString("F3", CultureInfo.InvariantCulture);
+
+        if (!hasPrevious)
+        {
+            WriteTestTime(filePath, elapsedMilliseconds);
+            Console.WriteLine($"FIRST TIME: {elapsedText} ms");
+            return;
+        }
+
+        string previousText = previousMilliseconds.ToString("F3", CultureInfo.InvariantCulture);
+        if (elapsedMilliseconds < previousMilliseconds)
+        {
+            WriteTestTime(filePath, elapsedMilliseconds);
+            Console.WriteLine($"BETTER TIME: {elapsedText} ms (previous {previousText} ms)");
+        }
+        else
+        {
+            Console.WriteLine($"WORSE TIME: {elapsedText} ms (best {previousText} ms)");
+        }
+    }
+
+    private static string GetTestTimeFilePath()
+    {
+        if (string.IsNullOrEmpty(_resultsDir))
+        {
+            return TestTimeFileName;
+        }
+
+        return Path.Combine(_resultsDir!, TestTimeFileName);
+    }
+
+    private static bool TryReadTestTime(string filePath, out double milliseconds)
+    {
+        if (!File.Exists(filePath))
+        {
+            milliseconds = 0.0;
+            return false;
+        }
+
+        string content = File.ReadAllText(filePath).Trim();
+        if (content.Length == 0)
+        {
+            milliseconds = 0.0;
+            return false;
+        }
+
+        return double.TryParse(content, NumberStyles.Float, CultureInfo.InvariantCulture, out milliseconds);
+    }
+
+    private static void WriteTestTime(string filePath, double milliseconds)
+    {
+        string? directory = Path.GetDirectoryName(filePath);
+        if (!string.IsNullOrEmpty(directory) && !Directory.Exists(directory))
+        {
+            Directory.CreateDirectory(directory);
+        }
+
+        string formatted = milliseconds.ToString("F6", CultureInfo.InvariantCulture);
+        File.WriteAllText(filePath, formatted);
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     internal static ulong CountOnes(ulong value)
     {
@@ -1419,6 +1618,8 @@ internal static class Program
                 _lastCompositeP = true;
                 return false;
             }
+
+            OnPrimeCandidateConfirmed();
         }
 
         searchedMersenne = true;

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -1,10 +1,12 @@
 using System.Buffers;
+using System.Collections.Generic;
 using System.Buffers.Binary;
 using System.Diagnostics;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
 using PerfectNumbers.Core.Gpu;
+using System.Collections.Concurrent;
 
 namespace PerfectNumbers.Core;
 
@@ -46,15 +48,22 @@ internal static partial class PrimeOrderCalculator
 
 	[ThreadStatic]
 	private static bool s_debugLoggingEnabled;
+	[ThreadStatic]
+	private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers;
 
 	[ThreadStatic]
-	private static List<ulong>? s_partialFactorPending;
+	private static ulong[]? s_smallPrimeStepperPreviousValues;
 
 	[ThreadStatic]
-	private static Stack<ulong>? s_partialFactorCompositeStack;
+	private static ulong[]? s_smallPrimeStepperRemainders;
 
 	[ThreadStatic]
-	private static Dictionary<ulong, int>? s_partialFactorCounts;
+	private static bool[]? s_smallPrimeStepperHasState;
+
+
+
+
+
 
 	private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
 
@@ -402,12 +411,13 @@ internal static partial class PrimeOrderCalculator
 			}
 
 			int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
-			List<ulong> candidates = new(capacity);
+			List<ulong> candidates = ThreadStaticPools.RentUlongList(capacity);
 			FactorEntry[] factorArray = orderFactors.Factors!;
 			// DebugLog("Building candidates list");
 			BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
 			if (candidates.Count == 0)
 			{
+				ThreadStaticPools.ReturnUlongList(candidates);
 				return false;
 			}
 
@@ -508,6 +518,7 @@ internal static partial class PrimeOrderCalculator
 						pool.Return(gpuPool, clearArray: false);
 					}
 
+					ThreadStaticPools.ReturnUlongList(candidates);
 					result = candidate;
 					return true;
 				}
@@ -520,6 +531,7 @@ internal static partial class PrimeOrderCalculator
 				index += batchSize;
 			}
 
+			ThreadStaticPools.ReturnUlongList(candidates);
 			// DebugLog("No candidate confirmed");
 			return false;
 		}
@@ -774,67 +786,32 @@ internal static partial class PrimeOrderCalculator
 		return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
 	}
 
-	private static List<ulong> AcquirePendingCompositeList(int capacityHint)
+	// [ThreadStatic]
-	{
+	private static ConcurrentDictionary<ulong, bool>? _primalityCache = [];
-		List<ulong>? list = s_partialFactorPending;
-		if (list is null)
-		{
-			list = new List<ulong>(capacityHint);
-			s_partialFactorPending = list;
-			return list;
-		}
 
-		list.Clear();
+	// [ThreadStatic]
-		if (list.Capacity < capacityHint)
+	// private static UInt128 _isPrimeHits;
-		{
-			list.Capacity = capacityHint;
-		}
-
-		return list;
-	}
-
-	private static Stack<ulong> AcquireCompositeStack(int capacityHint)
-	{
-		Stack<ulong>? stack = s_partialFactorCompositeStack;
-		if (stack is null)
-		{
-			stack = new Stack<ulong>(capacityHint);
-			s_partialFactorCompositeStack = stack;
-			return stack;
-		}
-
-		stack.Clear();
-		stack.EnsureCapacity(capacityHint);
-		return stack;
-	}
-
-	private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
-	{
-		Dictionary<ulong, int>? dictionary = s_partialFactorCounts;
-		if (dictionary is null)
-		{
-			dictionary = new Dictionary<ulong, int>(capacityHint);
-			s_partialFactorCounts = dictionary;
-			return dictionary;
-		}
-
-		dictionary.Clear();
-		dictionary.EnsureCapacity(capacityHint);
-		return dictionary;
-	}
 
 	private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config)
 	{
+		var primalityCache = _primalityCache;// ??= [];
 		if (value <= 1UL)
 		{
 			return PartialFactorResult.Empty;
 		}
 
 		const int FactorSlotCount = GpuSmallPrimeFactorSlots;
+
+		// stackalloc is faster than pooling
+		// ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount);
+		// int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount);
+		// Span<ulong> primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount);
+		// Span<int> exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount);
 		Span<ulong> primeSlots = stackalloc ulong[FactorSlotCount];
 		Span<int> exponentSlots = stackalloc int[FactorSlotCount];
-		primeSlots.Clear();
+		// We don't need to worry about leftovers, because we always use indexes within the calculated counts
-		exponentSlots.Clear();
+		// primeSlots.Clear();
+		// exponentSlots.Clear();
 
 		int factorCount = 0;
 		Dictionary<ulong, int>? counts = null;
@@ -849,190 +826,288 @@ internal static partial class PrimeOrderCalculator
 			gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
 		}
 
-		List<ulong> pending = AcquirePendingCompositeList(2);
+		List<ulong> pending = ThreadStaticPools.RentUlongList(2);
 		Stack<ulong>? compositeStack = null;
+		PartialFactorResult result;
 
-		try
+		if (!gpuFactored)
 		{
-			if (!gpuFactored)
+			// We don't need to worry about leftovers, because we always use indexes within the calculated counts
+			// primeSlots.Clear();
+			// exponentSlots.Clear();
+			if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
 			{
-				primeSlots.Clear();
+				useDictionary = true;
-				exponentSlots.Clear();
+				counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(factorCount, 8));
-				if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
+				remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
+			}
+		}
+
+		if (remaining > 1UL)
+		{
+			pending.Add(remaining);
+		}
+
+		if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
+		{
+			long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
+			compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count << 1, 4));
+			compositeStack.Push(remaining);
+			pending.Clear();
+
+			long timestamp = 0L;
+			while (compositeStack.Count > 0)
+			{
+				ulong composite = compositeStack.Pop();
+				// composite will never be smaller than 1 on the execution path
+				// if (composite <= 1UL)
+				// {
+				// 	continue;
+				// }
+
+				if (!primalityCache!.TryGetValue(composite, out bool isPrime))
 				{
-					useDictionary = true;
+					// Console.WriteLine($"Checking composite {composite}");				
-					counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8));
+					isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
-					remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
+					_ = primalityCache.TryAdd(composite, isPrime);
 				}
-			}
+				// else
+				// {
+				// 	Console.WriteLine($"isPrimeHits = {++_isPrimeHits}");
+				// }
 				
-			if (remaining > 1UL)
-			{
-				pending.Add(remaining);
-			}
-
-			if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
-			{
-				long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
-				compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4));
-				compositeStack.Push(remaining);
-				pending.Clear();
-
-				long timestamp = 0L; // reused for deadline checks.
-				while (compositeStack.Count > 0)
-				{
-					ulong composite = compositeStack.Pop();
-					if (composite <= 1UL)
-					{
-						continue;
-					}
-
-					bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
-					if (isPrime)
-					{
-						AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
-						continue;
-					}
-
-					timestamp = Stopwatch.GetTimestamp();
-					if (timestamp > deadlineTimestamp)
-					{
-						pending.Add(composite);
-						continue;
-					}
-
-					if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor))
-					{
-						pending.Add(composite);
-						continue;
-					}
-
-					ulong quotient = composite / factor;
-					compositeStack.Push(factor);
-					compositeStack.Push(quotient);
-				}
-			}
-
-			ulong cofactor = 1UL;
-			int pendingCount = pending.Count;
-			for (int i = 0; i < pendingCount; i++)
-			{
-				ulong composite = pending[i];
-				bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 				if (isPrime)
 				{
 					AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
-				}
-				else if (composite > 1UL)
-				{
-					cofactor = checked(cofactor * composite);
-				}
-			}
-
-			ArrayPool<FactorEntry> pool = ThreadStaticPools.FactorEntryPool;
-			if (useDictionary)
-			{
-				if ((counts is null || counts.Count == 0) && cofactor == value)
-				{
-					return PartialFactorResult.Rent(null, value, false, 0);
-				}
-
-				if (counts is null || counts.Count == 0)
-				{
-					return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
-				}
-
-				FactorEntry[] factors = pool.Rent(counts.Count);
-				int index = 0;
-				foreach (KeyValuePair<ulong, int> entry in counts)
-				{
-					factors[index] = new FactorEntry(entry.Key, entry.Value);
-					index++;
-				}
-
-				Span<FactorEntry> span = factors.AsSpan(0, index);
-				span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
-
-				bool fullyFactored = cofactor == 1UL;
-				return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index);
-			}
-
-			int actualCount = 0;
-			for (int i = 0; i < factorCount; i++)
-			{
-				if (primeSlots[i] != 0UL && exponentSlots[i] != 0)
-				{
-					actualCount++;
-				}
-			}
-
-			if (actualCount == 0)
-			{
-				if (cofactor == value)
-				{
-					return PartialFactorResult.Rent(null, value, false, 0);
-				}
-
-				return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
-			}
-
-			FactorEntry[] array = pool.Rent(actualCount);
-			int arrayIndex = 0;
-			for (int i = 0; i < factorCount; i++)
-			{
-				ulong primeValue = primeSlots[i];
-				int exponentValue = exponentSlots[i];
-				if (primeValue == 0UL || exponentValue == 0)
-				{
 					continue;
 				}
 
-				array[arrayIndex] = new FactorEntry(primeValue, exponentValue);
+				timestamp = Stopwatch.GetTimestamp();
-				arrayIndex++;
+				if (timestamp > deadlineTimestamp)
+				{
+					pending.Add(composite);
+					continue;
+				}
+
+				if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor))
+				{
+					pending.Add(composite);
+					continue;
+				}
+
+				ulong quotient = composite / factor;
+				compositeStack.Push(factor);
+				compositeStack.Push(quotient);
+			}
+		}
+
+		ulong cofactor = 1UL;
+		int pendingCount = pending.Count;
+		for (int i = 0; i < pendingCount; i++)
+		{
+			ulong composite = pending[i];
+			if (!primalityCache!.TryGetValue(composite, out bool isPrime))
+			{
+				// Console.WriteLine($"Checking composite {composite}");				
+				isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
+				_ = primalityCache.TryAdd(composite, isPrime);
+			}
+			if (isPrime)
+			{
+				AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
+			}
+			// composite is never smaller on the execution path
+			// else if (composite > 1UL)
+			else
+			{
+				cofactor = checked(cofactor * composite);
+			}
+		}
+
+		ArrayPool<FactorEntry> pool = ThreadStaticPools.FactorEntryPool;
+		if (useDictionary)
+		{
+			if ((counts is null || counts.Count == 0) && cofactor == value)
+			{
+				result = PartialFactorResult.Rent(null, value, false, 0);
+				goto ReturnResult;
 			}
 
-			Span<FactorEntry> arraySpan = array.AsSpan(0, arrayIndex);
+			if (counts is null || counts.Count == 0)
-			arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+			{
-			bool fullyFactoredArray = cofactor == 1UL;
+				result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
-			return PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex);
+				goto ReturnResult;
+			}
+
+			FactorEntry[] factors = pool.Rent(counts.Count);
+			int index = 0;
+			foreach (KeyValuePair<ulong, int> entry in counts)
+			{
+				factors[index] = new FactorEntry(entry.Key, entry.Value);
+				index++;
+			}
+
+			Span<FactorEntry> span = factors.AsSpan(0, index);
+			span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+
+			bool fullyFactored = cofactor == 1UL;
+			result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index);
+			goto ReturnResult;
 		}
-		finally
+
+		int actualCount = 0;
+		for (int i = 0; i < factorCount; i++)
 		{
-			pending.Clear();
+			if (primeSlots[i] != 0UL && exponentSlots[i] != 0)
-			counts?.Clear();
+			{
-			compositeStack?.Clear();
+				actualCount++;
+			}
 		}
+
+		if (actualCount == 0)
+		{
+			if (cofactor == value)
+			{
+				result = PartialFactorResult.Rent(null, value, false, 0);
+				goto ReturnResult;
+			}
+
+			result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
+			goto ReturnResult;
+		}
+
+		FactorEntry[] array = pool.Rent(actualCount);
+		int arrayIndex = 0;
+		for (int i = 0; i < factorCount; i++)
+		{
+			ulong primeValue = primeSlots[i];
+			int exponentValue = exponentSlots[i];
+			if (primeValue == 0UL || exponentValue == 0)
+			{
+				continue;
+			}
+
+			array[arrayIndex] = new FactorEntry(primeValue, exponentValue);
+			arrayIndex++;
+		}
+
+		Span<FactorEntry> arraySpan = array.AsSpan(0, arrayIndex);
+		arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+		bool fullyFactoredArray = cofactor == 1UL;
+		result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex);
+
+	ReturnResult:
+		// ThreadStaticPools.UlongPool.Return(primeSlotsArray);
+		// ThreadStaticPools.IntPool.Return(exponentSlotsArray);
+
+		if (counts is not null)
+		{
+			ThreadStaticPools.ReturnUlongIntDictionary(counts);
+		}
+
+		if (compositeStack is not null)
+		{
+			compositeStack.Clear();
+			ThreadStaticPools.ReturnUlongStack(compositeStack);
+		}
+
+		pending.Clear();
+		ThreadStaticPools.ReturnUlongList(pending);
+		return result;
+	}
+
+	private static CycleRemainderStepper[] EnsureSmallPrimeCycleSteppers(uint[] primes)
+	{
+		int primeCount = primes.Length;
+		CycleRemainderStepper[]? steppers = s_smallPrimeCycleSteppers;
+		if (steppers is null || steppers.Length < primeCount)
+		{
+			int length = primeCount;
+			steppers = new CycleRemainderStepper[length];
+			for (int i = 0; i < length; i++)
+			{
+				steppers[i] = new CycleRemainderStepper(primes[i]);
+			}
+
+			s_smallPrimeCycleSteppers = steppers;
+			s_smallPrimeStepperPreviousValues = new ulong[length];
+			s_smallPrimeStepperRemainders = new ulong[length];
+			s_smallPrimeStepperHasState = new bool[length];
+			return steppers;
+		}
+
+		if (s_smallPrimeStepperPreviousValues is null || s_smallPrimeStepperPreviousValues.Length != steppers.Length)
+		{
+			var previousValues = new ulong[steppers.Length];
+			if (s_smallPrimeStepperPreviousValues is { } existingPrevious)
+			{
+				System.Array.Copy(existingPrevious, previousValues, Math.Min(existingPrevious.Length, previousValues.Length));
+			}
+
+			s_smallPrimeStepperPreviousValues = previousValues;
+		}
+
+		if (s_smallPrimeStepperRemainders is null || s_smallPrimeStepperRemainders.Length != steppers.Length)
+		{
+			var remainders = new ulong[steppers.Length];
+			if (s_smallPrimeStepperRemainders is { } existingRemainders)
+			{
+				System.Array.Copy(existingRemainders, remainders, Math.Min(existingRemainders.Length, remainders.Length));
+			}
+
+			s_smallPrimeStepperRemainders = remainders;
+		}
+
+		if (s_smallPrimeStepperHasState is null || s_smallPrimeStepperHasState.Length != steppers.Length)
+		{
+			var hasState = new bool[steppers.Length];
+			if (s_smallPrimeStepperHasState is { } existingHasState)
+			{
+				System.Array.Copy(existingHasState, hasState, Math.Min(existingHasState.Length, hasState.Length));
+			}
+
+			s_smallPrimeStepperHasState = hasState;
+		}
+
+		return steppers;
 	}
 
 	private const int GpuSmallPrimeFactorSlots = 64;
 
 	private static bool TryPopulateSmallPrimeFactorsCpu(
-		ulong value,
+			ulong value,
-		uint limit,
+			uint limit,
-		Span<ulong> primeTargets,
+			Span<ulong> primeTargets,
-		Span<int> exponentTargets,
+			Span<int> exponentTargets,
-		out int factorCount,
+			out int factorCount,
-		out ulong remaining)
+			out ulong remaining)
 	{
 		factorCount = 0;
 		remaining = value;
 
-		if (value <= 1UL)
+		// value will never <= 1 in production code
-		{
+		// if (value <= 1UL)
-			return true;
+		// {
-		}
+		// 	return true;
+		// }
 
 		int capacity = Math.Min(primeTargets.Length, exponentTargets.Length);
-		if (capacity == 0)
+		// capacity will never equal 0 in production code
-		{
+		// if (capacity == 0)
-			return false;
+		// {
-		}
+		// 	return false;
+		// }
 
 		uint[] primes = PrimesGenerator.SmallPrimes;
 		ulong[] squares = PrimesGenerator.SmallPrimesPow2;
 		int primeCount = primes.Length;
 		ulong remainingLocal = value;
 		uint effectiveLimit = limit == 0 ? uint.MaxValue : limit;
+		CycleRemainderStepper[] cycleSteppers = EnsureSmallPrimeCycleSteppers(primes);
+		ulong[] stepperPreviousValues = s_smallPrimeStepperPreviousValues!;
+		ulong[] stepperRemainders = s_smallPrimeStepperRemainders!;
+		bool[] stepperHasState = s_smallPrimeStepperHasState!;
+		ulong valueForSteppers = value;
 
 		for (int i = 0; i < primeCount && remainingLocal > 1UL; i++)
 		{
@@ -1043,31 +1118,69 @@ internal static partial class PrimeOrderCalculator
 			}
 
 			ulong primeSquare = squares[i];
-			if (primeSquare != 0UL && primeSquare > remainingLocal)
+			// primeSquare will never == 0 in production code
+			// if (primeSquare != 0 && primeSquare > remainingLocal)
+			if (primeSquare > remainingLocal)
 			{
 				break;
 			}
 
-			ulong primeValue = primeCandidate;
+			// primeCandidate will never equal 0 in production code
-			if ((remainingLocal % primeValue) != 0UL)
+			// if (primeCandidate == 0UL)
+			// {
+			// 	continue;
+			// }
+
+			ref CycleRemainderStepper stepper = ref cycleSteppers[i];
+			ulong remainder;
+
+			if (!stepperHasState[i])
+			{
+				stepper.Reset();
+				remainder = stepper.Initialize(valueForSteppers);
+				stepperHasState[i] = true;
+			}
+			else
+			{
+				ulong previousValue = stepperPreviousValues[i];
+				if (valueForSteppers > previousValue)
+				{
+					remainder = stepper.ComputeNext(valueForSteppers);
+				}
+				else if (valueForSteppers < previousValue)
+				{
+					stepper.Reset();
+					remainder = stepper.Initialize(valueForSteppers);
+				}
+				else
+				{
+					remainder = stepperRemainders[i];
+				}
+			}
+
+			stepperPreviousValues[i] = valueForSteppers;
+			stepperRemainders[i] = remainder;
+
+			if (remainder != 0UL)
 			{
 				continue;
 			}
 
+			// TODO: Can we calculate this using residue stepping? We can try reusing modulo calculation for division result or the other way around
 			int exponent = 0;
 			do
 			{
-				remainingLocal /= primeValue;
+				remainingLocal /= primeCandidate;
 				exponent++;
 			}
-			while ((remainingLocal % primeValue) == 0UL);
+			while ((remainingLocal % primeCandidate) == 0UL);
 
 			if (factorCount >= capacity)
 			{
-				return false;
+				throw new InvalidOperationException($"Capacity is smaller than factor count");
 			}
 
-			primeTargets[factorCount] = primeValue;
+			primeTargets[factorCount] = primeCandidate;
 			exponentTargets[factorCount] = exponent;
 			factorCount++;
 		}
@@ -1118,14 +1231,18 @@ internal static partial class PrimeOrderCalculator
 
 	private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out ulong remaining)
 	{
-		Span<ulong> primeBuffer = stackalloc ulong[GpuSmallPrimeFactorSlots];
+		ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots);
-		Span<int> exponentBuffer = stackalloc int[GpuSmallPrimeFactorSlots];
+		int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots);
+		Span<ulong> primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
+		Span<int> exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
 		primeBuffer.Clear();
 		exponentBuffer.Clear();
 
 		if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _))
 		{
 			remaining = value;
+			ThreadStaticPools.UlongPool.Return(primeBufferArray);
+			ThreadStaticPools.IntPool.Return(exponentBufferArray);
 			return false;
 		}
 
@@ -1142,79 +1259,81 @@ internal static partial class PrimeOrderCalculator
 			counts[primeValue] = exponent;
 		}
 
+		ThreadStaticPools.UlongPool.Return(primeBufferArray);
+		ThreadStaticPools.IntPool.Return(exponentBufferArray);
 		return true;
 	}
-        private static long CreateDeadlineTimestamp(int milliseconds)
+	private static long CreateDeadlineTimestamp(int milliseconds)
-        {
+	{
-                if (milliseconds <= 0)
+		if (milliseconds <= 0)
-                {
+		{
-                        return long.MaxValue;
+			return long.MaxValue;
-                }
+		}
 
-                long stopwatchTicks = ConvertMillisecondsToStopwatchTicks(milliseconds);
+		long stopwatchTicks = ConvertMillisecondsToStopwatchTicks(milliseconds);
-                long startTimestamp = Stopwatch.GetTimestamp();
+		long startTimestamp = Stopwatch.GetTimestamp();
-                long deadline = startTimestamp + stopwatchTicks;
+		long deadline = startTimestamp + stopwatchTicks;
-                if (deadline < startTimestamp)
+		if (deadline < startTimestamp)
-                {
+		{
-                        return long.MaxValue;
+			return long.MaxValue;
-                }
+		}
 
-                return deadline;
+		return deadline;
-        }
+	}
 
-        private static long ConvertMillisecondsToStopwatchTicks(int milliseconds)
+	private static long ConvertMillisecondsToStopwatchTicks(int milliseconds)
-        {
+	{
-                long frequency = Stopwatch.Frequency;
+		long frequency = Stopwatch.Frequency;
-                long baseTicks = (frequency / 1000L) * milliseconds;
+		long baseTicks = (frequency / 1000L) * milliseconds;
-                long remainder = (frequency % 1000L) * milliseconds;
+		long remainder = (frequency % 1000L) * milliseconds;
-                return baseTicks + (remainder / 1000L);
+		return baseTicks + (remainder / 1000L);
-        }
+	}
 
-        private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor)
+	private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor)
-        {
+	{
-                factor = 0UL;
+		factor = 0UL;
-                if ((n & 1UL) == 0UL)
+		if ((n & 1UL) == 0UL)
-                {
+		{
-                        factor = 2UL;
+			factor = 2UL;
-                        return true;
+			return true;
-                }
+		}
 
-                long timestamp = 0L; // reused for deadline checks.
+		long timestamp = 0L; // reused for deadline checks.
-                while (true)
+		while (true)
-                {
+		{
-                        timestamp = Stopwatch.GetTimestamp();
+			timestamp = Stopwatch.GetTimestamp();
-                        if (timestamp > deadlineTimestamp)
+			if (timestamp > deadlineTimestamp)
-                        {
+			{
-                                return false;
+				return false;
-                        }
+			}
 
-                        ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+			ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
-                        ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
+			ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
-                        ulong y = x;
+			ulong y = x;
-                        ulong d = 1UL;
+			ulong d = 1UL;
 
-                        while (d == 1UL)
+			while (d == 1UL)
-                        {
+			{
-                                timestamp = Stopwatch.GetTimestamp();
+				timestamp = Stopwatch.GetTimestamp();
-                                if (timestamp > deadlineTimestamp)
+				if (timestamp > deadlineTimestamp)
-                                {
+				{
-                                        return false;
+					return false;
-                                }
+				}
 
-                                x = AdvancePolynomial(x, c, n);
+				x = AdvancePolynomial(x, c, n);
-                                y = AdvancePolynomial(y, c, n);
+				y = AdvancePolynomial(y, c, n);
-                                y = AdvancePolynomial(y, c, n);
+				y = AdvancePolynomial(y, c, n);
-                                ulong diff = x > y ? x - y : y - x;
+				ulong diff = x > y ? x - y : y - x;
-                                d = BinaryGcd(diff, n);
+				d = BinaryGcd(diff, n);
-                        }
+			}
 
-                        if (d != n)
+			if (d != n)
-                        {
+			{
-                                factor = d;
+				factor = d;
-                                return true;
+				return true;
-                        }
+			}
-                }
+		}
-        }
+	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus)
@@ -1318,7 +1437,7 @@ internal static partial class PrimeOrderCalculator
 			return;
 		}
 
-		counts ??= AcquireFactorCountDictionary(Math.Max(count, 8));
+		counts ??= ThreadStaticPools.RentUlongIntDictionary(Math.Max(count, 8));
 		CopyFactorsToDictionary(primes, exponents, count, counts);
 		count = 0;
 		useDictionary = true;

PerfectNumbers.Core/ThreadStaticPools.cs

@@ -79,6 +79,142 @@ namespace PerfectNumbers.Core
             }
         }
 
+        [ThreadStatic]
+        private static ArrayPool<int>? _intPool;
+
+        public static ArrayPool<int> IntPool
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get
+            {
+                return _intPool ??= ArrayPool<int>.Create();
+            }
+        }
+
+
+        [ThreadStatic]
+        private static List<List<ulong>>? _ulongListPool;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static List<ulong> RentUlongList(int capacityHint)
+        {
+            List<List<ulong>>? pool = _ulongListPool;
+            if (pool is not null && pool.Count > 0)
+            {
+                int lastIndex = pool.Count - 1;
+                List<ulong> list = pool[lastIndex];
+                pool.RemoveAt(lastIndex);
+                if (list.Count > 0)
+                {
+                    list.Clear();
+                }
+
+                if (list.Capacity < capacityHint)
+                {
+                    list.Capacity = capacityHint;
+                }
+
+                return list;
+            }
+
+            return new List<ulong>(capacityHint);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void ReturnUlongList(List<ulong> list)
+        {
+            List<List<ulong>>? pool = _ulongListPool;
+            if (pool is null)
+            {
+                pool = new List<List<ulong>>(4);
+                _ulongListPool = pool;
+            }
+
+            pool.Add(list);
+        }
+
+        [ThreadStatic]
+        private static List<Stack<ulong>>? _ulongStackPool;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Stack<ulong> RentUlongStack(int capacityHint)
+        {
+            List<Stack<ulong>>? pool = _ulongStackPool;
+            if (pool is not null && pool.Count > 0)
+            {
+                int lastIndex = pool.Count - 1;
+                Stack<ulong> stack = pool[lastIndex];
+                pool.RemoveAt(lastIndex);
+                if (stack.Count > 0)
+                {
+                    stack.Clear();
+                }
+
+                if (stack.Count < capacityHint)
+                {
+                    stack.EnsureCapacity(capacityHint);
+                }
+
+                return stack;
+            }
+
+            return new Stack<ulong>(capacityHint);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void ReturnUlongStack(Stack<ulong> stack)
+        {
+            List<Stack<ulong>>? pool = _ulongStackPool;
+            if (pool is null)
+            {
+                pool = new List<Stack<ulong>>(4);
+                _ulongStackPool = pool;
+            }
+
+            pool.Add(stack);
+        }
+
+        [ThreadStatic]
+        private static List<Dictionary<ulong, int>>? _ulongIntDictionaryPool;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Dictionary<ulong, int> RentUlongIntDictionary(int capacityHint)
+        {
+            List<Dictionary<ulong, int>>? pool = _ulongIntDictionaryPool;
+            if (pool is not null && pool.Count > 0)
+            {
+                int lastIndex = pool.Count - 1;
+                Dictionary<ulong, int> dictionary = pool[lastIndex];
+                pool.RemoveAt(lastIndex);
+                if (dictionary.Count > 0)
+                {
+                    dictionary.Clear();
+                }
+
+                if (dictionary.Count < capacityHint)
+                {
+                    dictionary.EnsureCapacity(capacityHint);
+                }
+
+                return dictionary;
+            }
+
+            return new Dictionary<ulong, int>(capacityHint);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void ReturnUlongIntDictionary(Dictionary<ulong, int> dictionary)
+        {
+            List<Dictionary<ulong, int>>? pool = _ulongIntDictionaryPool;
+            if (pool is null)
+            {
+                pool = new List<Dictionary<ulong, int>>(4);
+                _ulongIntDictionaryPool = pool;
+            }
+
+            pool.Add(dictionary);
+        }
+
         [ThreadStatic]
         private static Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? _mersenneFactorCacheDictionary;
 

PerfectNumbers.Core/ULongExtensions.cs

@@ -1,5 +1,4 @@
 using System;
-using System.Buffers;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using ILGPU.Algorithms;
@@ -346,55 +345,55 @@ public static class ULongExtensions
 
         ulong result = divisor.MontgomeryOne;
         ulong nPrime = divisor.NPrime;
-        ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount);
+
+        ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount);
+        Span<ulong> oddPowers = oddPowersArray.AsSpan(0, oddPowerCount);
+        InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers);
 
         int index = bitLength - 1;
-		while (index >= 0)
+        while (index >= 0)
-		{
-			if (((exponent >> index) & 1UL) == 0UL)
-			{
-				result = result.MontgomeryMultiply(result, modulus, nPrime);
-				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.MontgomeryMultiply(result, modulus, nPrime);
-			}
-
-			ulong mask = (1UL << windowLength) - 1UL;
-			ulong windowValue = (exponent >> windowStart) & mask;
-			int tableIndex = (int)((windowValue - 1UL) >> 1);
-			ulong multiplier = oddPowers[tableIndex];
-			result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
-
-			index = windowStart - 1;
-		}
-		
-		if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold)
-		{
-			ArrayPool<ulong>.Shared.Return(oddPowers);
-		}
-
-        if (keepMontgomery)
         {
-            return result;
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result = result.MontgomeryMultiply(result, modulus, nPrime);
+                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.MontgomeryMultiply(result, modulus, nPrime);
+            }
+
+            ulong mask = (1UL << windowLength) - 1UL;
+            ulong windowValue = (exponent >> windowStart) & mask;
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            ulong multiplier = oddPowers[tableIndex];
+            result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
+
+            index = windowStart - 1;
         }
 
-        return result.MontgomeryMultiply(1UL, modulus, nPrime);
+        ThreadStaticPools.UlongPool.Return(oddPowersArray);
+
+        if (!keepMontgomery)
+        {
+            result = result.MontgomeryMultiply(1UL, modulus, nPrime);
+        }
+
+        return result;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -533,26 +532,22 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+    private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> oddPowers)
     {
-		ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold ? new ulong[oddPowerCount] : ArrayPool<ulong>.Shared.Rent(oddPowerCount);
-		
         oddPowers[0] = divisor.MontgomeryTwo;
-        if (oddPowerCount == 1)
+        if (oddPowers.Length == 1)
         {
-            return oddPowers;
+            return;
         }
 
-        ulong previous,
+        ulong previous;
-			  square = divisor.MontgomeryTwoSquared;
+        ulong square = divisor.MontgomeryTwoSquared;
 
-        for (int i = 1; i < oddPowerCount; i++)
+        for (int i = 1; i < oddPowers.Length; i++)
         {
             previous = oddPowers[i - 1];
             oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
         }
-
-        return oddPowers;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]