Modulo operators replaced with steppers, where possible.

PerfectNumbers.Core/BidirectionalCycleRemainderStepper.cs

@@ -0,0 +1,62 @@
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+public struct BidirectionalCycleRemainderStepper(ulong modulus)
+{
+	private readonly ulong _modulus = modulus;
+	private ulong _previous = 0UL;
+	private ulong _remainder = 0UL;
+
+	public readonly ulong Modulus
+	{
+		[MethodImpl(MethodImplOptions.AggressiveInlining)]
+		get
+		{
+			return _modulus;
+		}
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public void Initialize(ulong value)
+	{
+		_remainder = value.ReduceCycleRemainder(_modulus);
+		_previous = value;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public ulong Step(ulong value)
+	{
+		long delta = unchecked((long)value - (long)_previous);
+		_previous = value;
+		if (delta == 0)
+		{
+			return _remainder;
+		}
+
+		if (delta > 0)
+		{
+			ulong step = ((ulong)delta).ReduceCycleRemainder(_modulus);
+			ulong sum = _remainder + step;
+			if (sum >= _modulus)
+			{
+				sum -= _modulus;
+			}
+
+			_remainder = sum;
+			return _remainder;
+		}
+
+		ulong backward = ((ulong)-delta).ReduceCycleRemainder(_modulus);
+		if (backward > _remainder)
+		{
+			ulong diff = backward - _remainder;
+			ulong wrap = diff.ReduceCycleRemainder(_modulus);
+			_remainder = wrap == 0 ? 0UL : _modulus - wrap;
+			return _remainder;
+		}
+
+		_remainder -= backward;
+		return _remainder;
+	}
+}

PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs

@@ -18,7 +18,7 @@ public struct BigIntegerCycleRemainderStepper(BigInteger cycleLength)
     public BigInteger Initialize(BigInteger value)
     {
         _previousValue = value;
-        _currentRemainder = value % _cycleLength;
+        _currentRemainder = value.ReduceCycleRemainder(_cycleLength);
         return _currentRemainder;
     }
 
@@ -27,11 +27,7 @@ public struct BigIntegerCycleRemainderStepper(BigInteger cycleLength)
     {
         BigInteger delta = value - _previousValue;
         _previousValue = value;
-        _currentRemainder += delta;
-        if (_currentRemainder >= _cycleLength)
-        {
-            _currentRemainder %= _cycleLength;
-        }
+        _currentRemainder = (_currentRemainder + delta).ReduceCycleRemainder(_cycleLength);
 
         return _currentRemainder;
     }

PerfectNumbers.Core/BigIntegerExtensions.cs

@@ -0,0 +1,113 @@
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+/// <summary>
+/// BigInteger helpers including cycle-remainder reduction using the subtract ladder that benchmarked best.
+/// Shared stats/warm set stay here so ulong reductions can reuse them through ULongExtensions.
+/// </summary>
+public static partial class BigIntegerExtensions
+{
+    // private const int UnrolledSubtractCount = 6;
+    // private static readonly HashSet<ulong> WarmModuli;
+    // private static long _hits;
+    // private static long _misses;
+
+    // static BigIntegerExtensions()
+    // {
+    //     ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
+    //     WarmModuli = new HashSet<ulong>(snapshot.Length);
+    //     for (int i = 0; i < snapshot.Length; i++)
+    //     {
+    //         ulong cycle = snapshot[i];
+    //         if (cycle > 1UL)
+    //         {
+    //             WarmModuli.Add(cycle);
+    //         }
+    //     }
+    // }
+
+    // public static void DumpCycleRemainderStats()
+    // {
+    //     long hits = Interlocked.Read(ref _hits);
+    //     long misses = Interlocked.Read(ref _misses);
+    //     Console.WriteLine($"CycleRemainderReducerCache: hits={hits}, misses={misses}");
+    // }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static BigInteger ReduceCycleRemainder(this BigInteger value, BigInteger modulus)
+    {
+        // Track(modulus);
+        // DumpCycleRemainderStats();
+
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        value -= modulus;
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        value -= modulus;
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        value -= modulus;
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        value -= modulus;
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        value -= modulus;
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        value -= modulus;
+        if (value < modulus)
+        {
+            return value;
+        }
+
+        return value % modulus;
+    }
+
+    // [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    // internal static void Track(in BigInteger modulus)
+    // {
+    //     if (modulus <= ulong.MaxValue)
+    //     {
+    //         Track((ulong)modulus);
+    //     }
+    //     else
+    //     {
+    //         Interlocked.Increment(ref _misses);
+    //     }
+    // }
+
+    // [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    // internal static void Track(ulong modulus)
+    // {
+    //     if (WarmModuli.Contains(modulus))
+    //     {
+    //         Interlocked.Increment(ref _hits);
+    //     }
+    //     else
+    //     {
+    //         Interlocked.Increment(ref _misses);
+    //     }
+    // }
+}

PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs

@@ -55,7 +55,7 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera
 
         // Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
         // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
-        var exponentStepper = new ExponentRemainderStepperCpu(divisorData);
+        var exponentStepper = new ExponentRemainderStepperCpu(divisorData, divisorCycle);
 
         var cycleStepper = new CycleRemainderStepper(divisorCycle);
 

PerfectNumbers.Core/CycleRemainderStepper.cs

@@ -19,14 +19,7 @@ public struct CycleRemainderStepper(ulong cycleLength)
     public ulong Initialize(ulong prime)
     {
         ulong cycleLength = _cycleLength;
-        ulong remainder = prime;
-        if (remainder >= cycleLength)
-        {
-            // TODO: Swap this `%` for the shared divisor-cycle remainder helper so initialization reuses the cached
-            // cycle deltas benchmarked faster than on-the-fly modulo work, matching the
-            // MersenneDivisorCycleLengthGpuBenchmarks winner for both CPU and GPU call sites.
-            remainder %= cycleLength;
-        }
+        ulong remainder = prime.ReduceCycleRemainder(cycleLength);
 
         _previousPrime = prime;
         _currentRemainder = remainder;
@@ -43,9 +36,7 @@ public struct CycleRemainderStepper(ulong cycleLength)
         if (ulong.MaxValue - _currentRemainder < delta)
         {
             UInt128 extended = (UInt128)_currentRemainder + delta;
-            // TODO: Replace this `%` with the UInt128-aware cycle reducer so large deltas use the cached subtraction
-            // ladder instead of falling back to the slower modulo implementation.
-            _currentRemainder = (ulong)(extended % cycleLength);
+            _currentRemainder = ((ulong)(extended % cycleLength)).ReduceCycleRemainder(cycleLength);
             return _currentRemainder;
         }
 
@@ -56,10 +47,7 @@ public struct CycleRemainderStepper(ulong cycleLength)
             _currentRemainder -= cycleLength;
             if (_currentRemainder >= cycleLength)
             {
-                // TODO: Route this `%` through the shared divisor-cycle helper so repeated wrap-arounds avoid
-                // modulo operations and match the benchmarked fast path highlighted in
-                // MersenneDivisorCycleLengthGpuBenchmarks.
-                _currentRemainder %= cycleLength;
+                _currentRemainder = _currentRemainder.ReduceCycleRemainder(cycleLength);
             }
         }
         return _currentRemainder;

PerfectNumbers.Core/ExponentRemainderStepperCpu.cs

@@ -1,3 +1,4 @@
+using System.Collections.Concurrent;
 using System.Runtime.CompilerServices;
 
 namespace PerfectNumbers.Core;
@@ -6,18 +7,20 @@ namespace PerfectNumbers.Core;
 /// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path.
 /// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order.
 /// </summary>
-internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor)
+internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL)
 {
-private readonly MontgomeryDivisorData _divisor = divisor;
+    private readonly MontgomeryDivisorData _divisor = divisor;
     private readonly ulong _modulus = divisor.Modulus;
     private readonly ulong _nPrime = divisor.NPrime;
     private readonly ulong _montgomeryOne = divisor.MontgomeryOne;
+    private readonly ulong _cycleLength = cycleLength;
     private ulong PreviousExponent = 0UL;
     private ulong _currentMontgomery = divisor.MontgomeryOne;
+    private static readonly ConcurrentDictionary<(ulong Cycle, ulong Delta), ulong> _cycleDeltaCache = new();
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor)
-		=> _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne;
+    public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor, ulong cycleLength)
+        => _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne && _cycleLength == cycleLength;
 
 	/// <summary>
 	/// Clears the cached Montgomery residue so the stepper can be reused for another divisor.
@@ -61,14 +64,14 @@ private readonly MontgomeryDivisorData _divisor = divisor;
     /// </summary>
     /// <param name="exponent">Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent.</param>
     /// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public ulong ComputeNextCpu(ulong exponent)
     {
         ulong delta = exponent - PreviousExponent;
-        // TODO: Once divisor cycle lengths are mandatory, pull the delta multiplier from the
-        // single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the
-        // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
-        ulong multiplier = delta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
+        ulong reducedDelta = ReduceDelta(delta);
+        ulong multiplier = reducedDelta == 0UL
+            ? _montgomeryOne
+            : reducedDelta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
         PreviousExponent = exponent;
         return ReduceCurrent();
@@ -85,9 +88,10 @@ private readonly MontgomeryDivisorData _divisor = divisor;
     public bool ComputeNextIsUnity(ulong exponent)
     {
         ulong delta = exponent - PreviousExponent;
-        // TODO: Reuse the divisor-cycle derived Montgomery delta once the cache exposes single-cycle
-        // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor.
-        ulong multiplier = delta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
+        ulong reducedDelta = ReduceDelta(delta);
+        ulong multiplier = reducedDelta == 0UL
+            ? _montgomeryOne
+            : reducedDelta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
         PreviousExponent = exponent;
         return _currentMontgomery == _montgomeryOne;
@@ -119,24 +123,62 @@ private readonly MontgomeryDivisorData _divisor = divisor;
     /// <param name="isUnity">Outputs whether the resulting residue equals Montgomery one.</param>
     /// <returns><see langword="true"/> when the advance succeeds.</returns>
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity)
-    {
-        // TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead
-        // of relying on the caller-provided delta so incremental scans remain in sync with the snapshot
-        // without computing additional cycles or mutating cache state.
-        _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(montgomeryDelta, _modulus, _nPrime);
+	public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity)
+	{
+        ulong delta = exponent - PreviousExponent;
+        ulong reducedDelta = ReduceDelta(delta);
+
+        // Prefer a cycle-reduced powmod so stepping stays aligned with the divisor-cycle cache.
+        ulong multiplier;
+        if (reducedDelta == 0UL)
+        {
+            multiplier = _montgomeryOne;
+        }
+        else if (_cycleLength != 0UL)
+        {
+            multiplier = reducedDelta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
+        }
+        else
+        {
+            multiplier = montgomeryDelta;
+        }
+
+        _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
         PreviousExponent = exponent;
         isUnity = _currentMontgomery == _montgomeryOne;
         return true;
     }
 
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private void InitializeCpuState(ulong exponent)
     {
-        _currentMontgomery = exponent.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
+        ulong rotation = ReduceDelta(exponent);
+        _currentMontgomery = rotation == 0UL
+            ? _montgomeryOne
+            : rotation.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
         PreviousExponent = exponent;
     }
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
     private readonly ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private readonly ulong ReduceDelta(ulong delta)
+    {
+        ulong cycleLength = _cycleLength;
+        if (cycleLength == 0UL)
+        {
+            return delta;
+        }
+
+        var key = (cycleLength, delta);
+        if (_cycleDeltaCache.TryGetValue(key, out ulong cached))
+        {
+            return cached;
+        }
+
+        ulong reduced = delta.ReduceCycleRemainder(cycleLength);
+        _cycleDeltaCache.TryAdd(key, reduced);
+        return reduced;
+    }
 }

PerfectNumbers.Core/HeuristicPrimeTester.cs

@@ -30,6 +30,23 @@ public sealed class HeuristicPrimeTester
 
 	private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
 
+	[ThreadStatic]
+	private static BidirectionalCycleRemainderStepper _mod10Stepper;
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static byte GetMod10Remainder(ulong value)
+	{
+		if (_mod10Stepper.Modulus != 10UL)
+		{
+			_mod10Stepper = new BidirectionalCycleRemainderStepper(10UL);
+			_mod10Stepper.Initialize(value);
+			return (byte)_mod10Stepper.Step(value);
+		}
+
+		// Caller must initialize once per thread before first use.
+		return (byte)_mod10Stepper.Step(value);
+	}
+
 	internal enum HeuristicDivisorGroup : byte
 	{
 		None = 0,
@@ -70,7 +87,7 @@ public sealed class HeuristicPrimeTester
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static bool IsPrimeCpu(ulong n)
 	{
-		byte nMod10 = (byte)n.Mod10();
+		byte nMod10 = GetMod10Remainder(n);
 		ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
 
 		if (maxDivisorSquare < 9UL)
@@ -85,7 +102,7 @@ public sealed class HeuristicPrimeTester
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n)
 	{
-		byte nMod10 = (byte)n.Mod10();
+		byte nMod10 = GetMod10Remainder(n);
 		ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
 		return IsPrimeGpu(gpu, n, maxDivisorSquare, nMod10);
 	}
@@ -132,11 +149,11 @@ public sealed class HeuristicPrimeTester
 						break;
 					}
 
-					ulong divisor = groupADivisors[groupAIndex];
-					if (n % divisor == 0UL)
-					{
-						return false;
-					}
+                    ulong divisor = groupADivisors[groupAIndex];
+                    if (n.ReduceCycleRemainder(divisor) == 0UL)
+                    {
+                        return false;
+                    }
 
 					groupAIndex++;
 					processed++;
@@ -186,10 +203,10 @@ public sealed class HeuristicPrimeTester
 					}
 
 					ulong divisor = groupADivisors[groupAIndex];
-					if (n % divisor == 0UL)
-					{
-						return false;
-					}
+                    if (n.ReduceCycleRemainder(divisor) == 0UL)
+                    {
+                        return false;
+                    }
 
 					groupAIndex++;
 					processed++;
@@ -220,7 +237,7 @@ public sealed class HeuristicPrimeTester
 						break;
 					}
 
-					if (n % divisor == 0UL)
+					if (n.ReduceCycleRemainder(divisor) == 0UL)
 					{
 						return false;
 					}
@@ -241,7 +258,6 @@ public sealed class HeuristicPrimeTester
 
 	private static bool HeuristicTrialDivisionGpuDetectsDivisor(PrimeOrderCalculatorAccelerator gpu, ulong n, ulong maxDivisorSquare, byte nMod10)
 	{
-		// GpuPrimeWorkLimiter.Acquire();
 		int acceleratorIndex = gpu.AcceleratorIndex;
 		var kernelLauncher = gpu.HeuristicCombinedTrialDivisionKernelLauncher;
 		bool compositeDetected;
@@ -252,13 +268,6 @@ public sealed class HeuristicPrimeTester
 		var stream = AcceleratorStreamPool.Rent(acceleratorIndex);
 		flagView1DView.CopyFromCPU(stream, ref compositeFlag, 1);
 
-		//         Index1D index,
-        // ArrayView<int> resultFlag,
-        // ulong n,
-        // ulong maxDivisorSquare,
-        // HeuristicGpuDivisorTableKind tableKind,
-        // HeuristicGpuDivisorTables tables)
-
 		kernelLauncher(
 				stream,
 				1,
@@ -271,48 +280,9 @@ public sealed class HeuristicPrimeTester
 		flagView1DView.CopyToCPU(stream, ref compositeFlag, 1);
 		stream.Synchronize();
 
-		compositeDetected = compositeFlag != 0;
-		// compositeDetected = compositeFlag != 0;
-
-		// if (!compositeDetected)
-		// {
-		// 	ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
-		// 	for (int i = 0; i < endingOrder.Length && !compositeDetected; i++)
-		// 	{
-		// 		byte ending = endingOrder[i];
-		// 		HeuristicGpuDivisorTableKind tableKind = ending switch
-		// 		{
-		// 			1 => HeuristicGpuDivisorTableKind.GroupBEnding1,
-		// 			7 => HeuristicGpuDivisorTableKind.GroupBEnding7,
-		// 			9 => HeuristicGpuDivisorTableKind.GroupBEnding9,
-		// 			_ => HeuristicGpuDivisorTableKind.GroupA,
-		// 		};
-
-		// 		if (tableKind == HeuristicGpuDivisorTableKind.GroupA)
-		// 		{
-		// 			continue;
-		// 		}
-
-		// 		int divisorLength = GetGroupBDivisors(ending).Length;
-
-		// 		compositeFlag = 0;
-		// 		flagView1D.CopyFromCPU(stream, ref compositeFlag, 1);
-		// 		kernelLauncher(
-		// 			stream,
-		// 				1,
-		// 				flagView1D,
-		// 				n,
-		// 				maxDivisorSquare,
-		// 				tableKind,
-		// 				gpu.HeuristicGpuTables);
-		// 		flagView1D.CopyToCPU(stream, ref compositeFlag, 1);
-		// 		stream.Synchronize();
-		// 		compositeDetected = compositeFlag != 0;
-		// 	}
-		// }
-
 		AcceleratorStreamPool.Return(acceleratorIndex, stream);
-		// GpuPrimeWorkLimiter.Release();
+		compositeDetected = compositeFlag != 0;
+
 		return compositeDetected;
 	}
 
@@ -556,6 +526,10 @@ public sealed class HeuristicPrimeTester
 		private readonly ReadOnlySpan<ulong> groupADivisorSquares;
 		private int groupAIndex;
 		private Span<HeuristicGroupBSequenceState> groupBStates;
+		private CycleRemainderStepper endingStepper;
+		private CycleRemainderStepper wheelStepper;
+		private bool hasRemainders;
+		private ulong lastDivisor;
 
 		public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
 		{
@@ -566,6 +540,10 @@ public sealed class HeuristicPrimeTester
 
 			int count = InitializeGroupBStates(nMod10, groupBBuffer);
 			groupBStates = count == 0 ? Span<HeuristicGroupBSequenceState>.Empty : groupBBuffer[..count];
+			endingStepper = new CycleRemainderStepper(10UL);
+			wheelStepper = new CycleRemainderStepper(Wheel210);
+			hasRemainders = false;
+			lastDivisor = 0UL;
 		}
 
 		public bool TryGetNext(out HeuristicDivisorCandidate candidate)
@@ -586,13 +564,14 @@ public sealed class HeuristicPrimeTester
 				HeuristicDivisorGroup group = currentIndex < GroupAConstantCount
 					? HeuristicDivisorGroup.GroupAConstant
 					: HeuristicDivisorGroup.GroupAWheel;
+				(byte ending, ushort wheelResidue) = UpdateRemainders(divisor);
 
 				candidate = new HeuristicDivisorCandidate(
 					divisor,
 					group,
-					(byte)(divisor % 10UL),
+					ending,
 					0,
-					(ushort)(divisor % Wheel210));
+					wheelResidue);
 				return true;
 			}
 
@@ -623,14 +602,35 @@ public sealed class HeuristicPrimeTester
 
 			ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex];
 			bestState.Advance();
+			(byte candidateEnding, ushort wheelResidueBest) = UpdateRemainders(bestCandidate);
 			candidate = new HeuristicDivisorCandidate(
 				bestCandidate,
 				HeuristicDivisorGroup.GroupB,
-				bestState.Ending,
+				candidateEnding,
 				bestState.PriorityIndex,
-				(ushort)(bestCandidate % Wheel210));
+				wheelResidueBest);
 			return true;
 		}
+
+		private (byte ending, ushort wheelResidue) UpdateRemainders(ulong divisor)
+		{
+			byte ending;
+			ushort wheelResidue;
+			if (!hasRemainders)
+			{
+				ending = (byte)endingStepper.Initialize(divisor);
+				wheelResidue = (ushort)wheelStepper.Initialize(divisor);
+				hasRemainders = true;
+			}
+			else
+			{
+				ending = (byte)endingStepper.ComputeNext(divisor);
+				wheelResidue = (ushort)wheelStepper.ComputeNext(divisor);
+			}
+
+			lastDivisor = divisor;
+			return (ending, wheelResidue);
+		}
 	}
 
 
@@ -684,6 +684,10 @@ public sealed class HeuristicPrimeTester
 		private readonly GpuUInt128 step;
 		private readonly GpuUInt128 limit;
 		private GpuUInt128 current;
+		private CycleRemainderStepper endingStepper;
+		private CycleRemainderStepper wheelStepper;
+		private byte currentEnding;
+		private ushort currentWheelResidue;
 		private MersenneDivisorResidueStepper residueStepper;
 		private bool active;
 		private ulong processedCount;
@@ -705,6 +709,21 @@ public sealed class HeuristicPrimeTester
 
 			current = hasCandidates ? firstDivisor : GpuUInt128.Zero;
 			residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default;
+			if (hasCandidates)
+			{
+				endingStepper = new CycleRemainderStepper(10UL);
+				currentEnding = (byte)endingStepper.Initialize(firstDivisor.Low);
+
+				wheelStepper = new CycleRemainderStepper(Wheel210);
+				currentWheelResidue = (ushort)wheelStepper.Initialize(firstDivisor.Low);
+			}
+			else
+			{
+				endingStepper = default;
+				wheelStepper = default;
+				currentEnding = 0;
+				currentWheelResidue = 0;
+			}
 			active = hasCandidates;
 			processedCount = 0UL;
 			lastDivisor = 0UL;
@@ -724,7 +743,7 @@ public sealed class HeuristicPrimeTester
 
 				if (admissible)
 				{
-					candidate = CreateCandidate(value);
+					candidate = CreateCandidate(value, currentEnding, currentWheelResidue);
 					return true;
 				}
 			}
@@ -746,6 +765,8 @@ public sealed class HeuristicPrimeTester
 
 			current = next;
 			residueStepper.Advance();
+			currentEnding = (byte)endingStepper.ComputeNext(next.Low);
+			currentWheelResidue = (ushort)wheelStepper.ComputeNext(next.Low);
 		}
 
 		public readonly ulong ProcessedCount => processedCount;
@@ -754,9 +775,8 @@ public sealed class HeuristicPrimeTester
 
 		public readonly bool Exhausted => !active;
 
-		private static HeuristicDivisorCandidate CreateCandidate(ulong value)
+		private static HeuristicDivisorCandidate CreateCandidate(ulong value, byte ending, ushort wheelResidue)
 		{
-			byte ending = (byte)(value % 10UL);
 			HeuristicDivisorGroup group;
 
 			if (value == 3UL || value == 7UL || value == 11UL || value == 13UL)
@@ -772,8 +792,7 @@ public sealed class HeuristicPrimeTester
 				group = HeuristicDivisorGroup.GroupB;
 			}
 
-			ushort residue = (ushort)(value % Wheel210);
-			return new HeuristicDivisorCandidate(value, group, ending, 0, residue);
+			return new HeuristicDivisorCandidate(value, group, ending, 0, wheelResidue);
 		}
 	}
 }

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -266,7 +266,7 @@ internal static partial class PrimeOrderCalculator
             int exponent = entries[i].Value;
             for (int iteration = 0; iteration < exponent; iteration++)
             {
-                if ((order % primeFactor) != UInt128.Zero)
+                if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero)
                 {
                     UInt128 candidate = order / primeFactor;
                     if (Pow2ModWideCpu(candidate, prime, divisorData) == UInt128.One)
@@ -339,14 +339,13 @@ internal static partial class PrimeOrderCalculator
 
         buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value));
 
-		// TODO: Investigating using stepper to replace modulo operation
         for (int i = 0; i < length; i++)
         {
             UInt128 primeFactor = buffer[i].Value;
             int exponent = buffer[i].Exponent;
             for (int iteration = 0; iteration < exponent; iteration++)
             {
-                if ((order % primeFactor) == UInt128.Zero)
+                if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero)
                 {
                     UInt128 reduced = order / primeFactor;
                     if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One)
@@ -405,7 +404,7 @@ internal static partial class PrimeOrderCalculator
             UInt128 reduced = order;
             for (int iteration = 0; iteration < span[i].Exponent; iteration++)
             {
-                if ((reduced % primeFactor) != UInt128.Zero)
+                if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero)
                 {
                     break;
                 }
@@ -2583,7 +2582,7 @@ internal static partial class PrimeOrderCalculator
 			// 	continue;
 			// }
 
-			if ((remainingLocal % primeValue) != 0UL)
+			if (remainingLocal.ReduceCycleRemainder(primeValue) != 0UL)
 			{
 				continue;
 			}
@@ -2745,16 +2744,16 @@ internal static partial class PrimeOrderCalculator
 	private static ulong AdvancePolynomialCpu(ulong x, ulong c, ulong modulus)
 	{
 		UInt128 value = (UInt128)x * x + c;
-		return (ulong)(value % modulus);
+		return (ulong)value.ReduceCycleRemainder(modulus);
 	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	private static ulong AdvancePolynomialTwiceCpu(ulong x, ulong c, ulong modulus)
 	{
 		UInt128 value = (UInt128)x * x + c;
-		var remainder = (ulong)(value % modulus);
+		var remainder = (ulong)value.ReduceCycleRemainder(modulus);
 		value = (UInt128)remainder * remainder + c;
-		return (ulong)(value % modulus);
+		return (ulong)value.ReduceCycleRemainder(modulus);
 	}
 
 	private static ulong BinaryGcd(ulong a, ulong b)
@@ -2897,7 +2896,7 @@ internal static partial class PrimeOrderCalculator
 			int exponent = entries[i].Value;
 			for (int iteration = 0; iteration < exponent; iteration++)
 			{
-				ulong remainder = order % primeFactor;
+				ulong remainder = order.ReduceCycleRemainder(primeFactor);
 				if (remainder != 0UL)
 				{
 					break;
@@ -2951,7 +2950,7 @@ internal static partial class PrimeOrderCalculator
 		{
 			int exponent = 1;
 			ulong remaining = quotient;
-			while ((remaining % factor) == 0UL)
+			while (remaining.ReduceCycleRemainder(factor) == 0UL)
 			{
 				remaining /= factor;
 				exponent++;
@@ -3062,7 +3061,7 @@ internal static partial class PrimeOrderCalculator
             UInt128 reduced = candidate;
             for (int iteration = 0; iteration < span[i].Exponent; iteration++)
             {
-                if ((reduced % primeFactor) != UInt128.Zero)
+                if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero)
                 {
                     break;
                 }

PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs

@@ -307,19 +307,19 @@ internal static partial class PrimeOrderCalculator
         {
             UInt128 primeFactor = span[i].Value;
             UInt128 reduced = order;
-            for (int iteration = 0; iteration < span[i].Exponent; iteration++)
-            {
-                if ((reduced % primeFactor) != UInt128.Zero)
-                {
-                    break;
-                }
+			for (int iteration = 0; iteration < span[i].Exponent; iteration++)
+			{
+				if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero)
+				{
+					break;
+				}
 
-                reduced /= primeFactor;
-                if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One)
-                {
-                    return false;
-                }
-            }
+				reduced /= primeFactor;
+				if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One)
+				{
+					return false;
+				}
+			}
         }
 
         return true;

PerfectNumbers.Core/PrimeTester.cs

@@ -1,7 +1,5 @@
-using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
-using ILGPU;
 using ILGPU.Runtime;
 using PerfectNumbers.Core.Gpu.Accelerators;
 
@@ -72,7 +70,7 @@ public sealed class PrimeTester
                 }
 
                 divisor = smallPrimeDivisors[i];
-                uint deltaMod = delta >= divisor ? (uint)(delta % divisor) : (uint)delta;
+                uint deltaMod = delta.ReduceCycleRemainder(divisor);
                 remainder = remainders[i] + deltaMod;
                 if (remainder >= divisor)
                 {
@@ -105,7 +103,7 @@ public sealed class PrimeTester
             }
 
             divisor = smallPrimeDivisors[i];
-            remainder = (uint)(n % divisor);
+            remainder = n.ReduceCycleRemainder(divisor);
             remainders[i] = remainder;
 
             if (remainder == 0UL)

PerfectNumbers.Core/PrimeTesterByLastDigit.cs

@@ -1,11 +1,29 @@
 using System.Runtime.CompilerServices;
-using ILGPU.Runtime;
-using PerfectNumbers.Core.Gpu.Accelerators;
 
 namespace PerfectNumbers.Core;
 
 public sealed class PrimeTesterByLastDigit
 {
+	[ThreadStatic]
+	private static BidirectionalCycleRemainderStepper _mod10Stepper;
+
+	[ThreadStatic]
+	private static bool _mod10Initialized;
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static byte GetMod10(ulong value)
+	{
+		if (!_mod10Initialized || _mod10Stepper.Modulus != 10UL)
+		{
+			_mod10Stepper = new BidirectionalCycleRemainderStepper(10UL);
+			_mod10Stepper.Initialize(value);
+			_mod10Initialized = true;
+			return (byte)_mod10Stepper.Step(value);
+		}
+
+		return (byte)_mod10Stepper.Step(value);
+	}
+
 	public static bool IsPrimeCpu(ulong n)
 	{
 		// The below IFs never trigger in production code of EvenPerfectBitScanner on --mersenne=bydivisor path.
@@ -27,7 +45,7 @@ public sealed class PrimeTesterByLastDigit
 
 		uint[] smallPrimeDivisors;
 		ulong[] smallPrimeDivisorsMul;
-		byte nMod10 = n.Mod10();
+		byte nMod10 = GetMod10(n);
 
 		switch (nMod10)
 		{

PerfectNumbers.Core/ThreadStaticPools.cs

@@ -360,9 +360,9 @@ namespace PerfectNumbers.Core
         private static bool _hasExponentRemainderStepperCpu;
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData)
+        internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData, ulong cycleLength = 0UL)
         {
-            if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData))
+            if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData, cycleLength))
             {
                 ExponentRemainderStepperCpu stepper = _exponentRemainderStepperCpu;
                 _hasExponentRemainderStepperCpu = false;
@@ -370,7 +370,7 @@ namespace PerfectNumbers.Core
                 return stepper;
             }
 
-            return new ExponentRemainderStepperCpu(divisorData);
+            return new ExponentRemainderStepperCpu(divisorData, cycleLength);
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]

PerfectNumbers.Core/UInt128Extensions.cs

@@ -306,6 +306,53 @@ public static class UInt128Extensions
 		return remainder.Pow2MontgomeryModWindowed(modulus);
 	}
 
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 ReduceCycleRemainder(this UInt128 value, UInt128 modulus)
+	{
+		if (modulus == UInt128.Zero || value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		return value % modulus;
+	}
+
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	private static bool ShouldUseSingleBit(UInt128 exponent) => (exponent >> 64) == UInt128.Zero && (ulong)exponent <= Pow2WindowFallbackThreshold;
 

PerfectNumbers.Core/UIntExtensions.cs

@@ -13,9 +13,59 @@ public static class UIntExtensions
     public static uint Mod3(uint value) => value - (uint)((value * Mod3Mask) >> 33) * 3U;
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static uint Mod5(this uint value) => value - (uint)((value * Mod5Mask) >> 34) * 5U;
+	public static uint Mod5(this uint value) => value - (uint)((value * Mod5Mask) >> 34) * 5U;
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static uint Mod8(this uint value) => value & 7U;
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static uint Mod8(this uint value) => value & 7U;
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static uint ReduceCycleRemainder(this uint value, uint modulus)
+	{
+		if (modulus == 0U || value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		return value % modulus;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static uint ReduceCycleRemainder(this ulong value, uint modulus)
+		=> (uint)value.ReduceCycleRemainder((ulong)modulus);
 }
-

PerfectNumbers.Core/ULongExtensions.cs

@@ -7,8 +7,6 @@ namespace PerfectNumbers.Core;
 public static partial class ULongExtensions
 {
 	private const ulong Pow2WindowFallbackThreshold = 32UL;
-	private const int Pow2WindowOddPowerCount = 1 << (PerfectNumberConstants.Pow2WindowSize - 1);
-
 	private static readonly UInt128 Pow2WindowedModulusThreshold = (UInt128)PerfectNumberConstants.MaxQForDivisorCycles;
 
 	private static readonly byte[] BitStatsZeroCount = new byte[256];
@@ -1135,6 +1133,56 @@ public static partial class ULongExtensions
 		return x1;
 	}
 
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong ReduceCycleRemainder(this ulong value, ulong modulus)
+	{
+		// BigIntegerExtensions.Track(modulus);
+		// BigIntegerExtensions.DumpCycleRemainderStats();
+
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		value -= modulus;
+		if (value < modulus)
+		{
+			return value;
+		}
+
+		return value % modulus;
+	}
+
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static bool SharesFactorWithExponentMinusOneCpu(this ulong exponent)
 	{