Merge pull request #92 from mlemanczyk/codex/refactor-cycleremainderstepper-and-mersennenumberdivisor

Initialize cycle remainder stepping upfront

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -307,7 +307,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 return;
             }
 
-            bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL;
+            bool useDivisorCycles = _owner._useDivisorCycles;
+            if (useDivisorCycles && divisorCycle == 0UL)
+            {
+                throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}.");
+            }
+
+            bool cycleEnabled = useDivisorCycles;
             MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
 
             var exponentStepper = new ExponentRemainderStepper(divisorData);
@@ -321,9 +327,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             {
                 var stepper = new CycleRemainderStepper(divisorCycle);
 
-                for (int i = 0; i < length; i++)
+                ulong remainder = stepper.Initialize(primes[0]);
+                hits[0] = remainder == 0UL
+                    ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0)
+                    : (byte)0;
+
+                for (int i = 1; i < length; i++)
                 {
-                    ulong remainder = stepper.ComputeNext(primes[i]);
+                    remainder = stepper.ComputeNext(primes[i]);
                     if (remainder != 0UL)
                     {
                         hits[i] = 0;

PerfectNumbers.Core/CycleRemainderStepper.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Numerics;
 
 namespace PerfectNumbers.Core;
@@ -11,14 +12,17 @@ internal struct CycleRemainderStepper
 
     public CycleRemainderStepper(ulong cycleLength)
     {
+        if (cycleLength == 0UL)
+        {
+            throw new ArgumentOutOfRangeException(nameof(cycleLength), "Cycle length must be non-zero.");
+        }
+
         _cycleLength = cycleLength;
         _previousPrime = 0UL;
         _currentRemainder = 0UL;
         _hasState = false;
     }
 
-    public bool HasCycle => _cycleLength != 0UL;
-
     public void Reset()
     {
         _previousPrime = 0UL;
@@ -26,52 +30,53 @@ internal struct CycleRemainderStepper
         _hasState = false;
     }
 
-    public void InitializeState(ulong previousPrime, ulong previousRemainder, bool hasState)
+    public ulong Initialize(ulong prime)
     {
-        _previousPrime = previousPrime;
+        ulong cycleLength = _cycleLength;
-        _currentRemainder = previousRemainder;
+        ulong remainder = prime;
-        _hasState = hasState;
+        if (remainder >= cycleLength)
+        {
+            remainder %= cycleLength;
+        }
+        _previousPrime = prime;
+        _currentRemainder = remainder;
+        _hasState = true;
+        return remainder;
     }
 
     public ulong ComputeNext(ulong prime)
     {
-        if (_cycleLength == 0UL)
+        if (!_hasState)
         {
-            _previousPrime = prime;
+            throw new InvalidOperationException("CycleRemainderStepper must be initialized before computing.");
-            _currentRemainder = 0UL;
-            _hasState = true;
-            return 0UL;
         }
 
-        if (!_hasState || prime <= _previousPrime)
+        if (prime <= _previousPrime)
         {
-            _currentRemainder = prime % _cycleLength;
+            throw new ArgumentOutOfRangeException(nameof(prime), "Primes must be processed in strictly increasing order.");
-            _hasState = true;
-        }
-        else
-        {
-            ulong delta = prime - _previousPrime;
-            if (ulong.MaxValue - _currentRemainder < delta)
-            {
-                UInt128 extended = (UInt128)_currentRemainder + delta;
-                _currentRemainder = (ulong)(extended % _cycleLength);
-            }
-            else
-            {
-                _currentRemainder += delta;
-
-                if (_currentRemainder >= _cycleLength)
-                {
-                    _currentRemainder -= _cycleLength;
-                    if (_currentRemainder >= _cycleLength)
-                    {
-                        _currentRemainder %= _cycleLength;
-                    }
-                }
-            }
         }
 
+        ulong cycleLength = _cycleLength;
+        ulong delta = prime - _previousPrime;
         _previousPrime = prime;
+
+        if (ulong.MaxValue - _currentRemainder < delta)
+        {
+            UInt128 extended = (UInt128)_currentRemainder + delta;
+            _currentRemainder = (ulong)(extended % cycleLength);
+            return _currentRemainder;
+        }
+
+        _currentRemainder += delta;
+
+        if (_currentRemainder >= cycleLength)
+        {
+            _currentRemainder -= cycleLength;
+            if (_currentRemainder >= cycleLength)
+            {
+                _currentRemainder %= cycleLength;
+            }
+        }
         return _currentRemainder;
     }
 

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -391,8 +391,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                                 return;
                         }
 
-                        EnsureCapacity(gpuBatchSize);
-
                         MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
                         var exponentStepper = new ExponentRemainderStepper(divisorData);
                         if (!exponentStepper.IsValidModulus)
@@ -407,7 +405,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                                 throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}.");
                         }
 
-                        CycleRemainderStepper cycleStepper = useDivisorCycles ? new CycleRemainderStepper(divisorCycle) : default;
+                        CycleRemainderStepper cycleStepper = default;
+                        ulong initialCycleRemainder = 0UL;
+                        bool cycleRemainderPending = false;
+                        if (useDivisorCycles)
+                        {
+                                cycleStepper = new CycleRemainderStepper(divisorCycle);
+                                initialCycleRemainder = cycleStepper.Initialize(primes[0]);
+                                cycleRemainderPending = true;
+                        }
 
                         Accelerator accelerator = _accelerator;
                         Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel;
@@ -437,7 +443,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
                                         if (useDivisorCycles)
                                         {
-                                                ulong remainder = cycleStepper.ComputeNext(primeValue);
+                                                ulong remainder;
+                                                if (cycleRemainderPending)
+                                                {
+                                                        remainder = initialCycleRemainder;
+                                                        cycleRemainderPending = false;
+                                                }
+                                                else
+                                                {
+                                                        remainder = cycleStepper.ComputeNext(primeValue);
+                                                }
                                                 if (remainder != 0UL)
                                                 {
                                                         continue;
@@ -507,43 +522,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         }
                 }
 
-                private void EnsureCapacity(int requiredCapacity)
+                public void Dispose()
                 {
-                        if (requiredCapacity <= _capacity)
+                        if (_disposed)
                         {
                                 return;
-                        }
-
-                        _exponentsBuffer.Dispose();
-                        _resultsBuffer.Dispose();
-                        ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
-                        ArrayPool<int>.Shared.Return(_positionBuffer, clearArray: false);
-                        ArrayPool<byte>.Shared.Return(_modeBuffer, clearArray: false);
-
-                        int newCapacity = _capacity;
-                        if (newCapacity < 1)
-                        {
-                                newCapacity = 1;
-                        }
-
-                        while (newCapacity < requiredCapacity)
-                        {
-                                newCapacity <<= 1;
-                        }
-
-                        _exponentsBuffer = _accelerator.Allocate1D<ulong>(newCapacity);
-                        _resultsBuffer = _accelerator.Allocate1D<ulong>(newCapacity);
-                        _hostBuffer = ArrayPool<ulong>.Shared.Rent(newCapacity);
-                        _positionBuffer = ArrayPool<int>.Shared.Rent(newCapacity);
-                        _modeBuffer = ArrayPool<byte>.Shared.Rent(newCapacity);
-                        _capacity = newCapacity;
-                }
-
-		public void Dispose()
-		{
-			if (_disposed)
-			{
-				return;
 			}
 
                         _disposed = true;