Trim GPU exponent batches using cached divisor cycles

EvenPerfectBitScanner/Program.cs

@@ -821,6 +821,8 @@ internal static class Program
                                 return;
                         }
 
+                        states.Sort(static (left, right) => left.Prime.CompareTo(right.Prime));
+
                         int stateCount = states.Count;
                         ulong[] primeValues = new ulong[stateCount];
                         ulong[] allowedMaxValues = new ulong[stateCount];

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -310,6 +310,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL;
             MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
 
+            var exponentStepper = new ExponentRemainderStepper(divisorData);
+            if (!exponentStepper.IsValidModulus)
+            {
+                hits.Clear();
+                return;
+            }
+
             if (cycleEnabled)
             {
                 var stepper = new CycleRemainderStepper(divisorCycle);
@@ -323,8 +330,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                         continue;
                     }
 
-                    ulong residue = remainder.Pow2MontgomeryModFromCycleRemainder(divisorData);
+                    hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0;
-                    hits[i] = residue == 1UL ? (byte)1 : (byte)0;
                 }
 
                 return;
@@ -332,8 +338,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
             for (int i = 0; i < length; i++)
             {
-                ulong residue = primes[i].Pow2MontgomeryMod(divisorData);
+                hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0;
-                hits[i] = residue == 1UL ? (byte)1 : (byte)0;
             }
         }
 

PerfectNumbers.Core/ExponentRemainderStepper.cs

@@ -0,0 +1,124 @@
+namespace PerfectNumbers.Core;
+
+internal struct ExponentRemainderStepper
+{
+    private readonly MontgomeryDivisorData _divisor;
+    private readonly ulong _modulus;
+    private readonly ulong _nPrime;
+    private readonly ulong _montgomeryOne;
+    private ulong _previousExponent;
+    private ulong _currentMontgomery;
+    private bool _hasState;
+
+    public ExponentRemainderStepper(in MontgomeryDivisorData divisor)
+    {
+        _divisor = divisor;
+        _modulus = divisor.Modulus;
+        _nPrime = divisor.NPrime;
+        _montgomeryOne = divisor.MontgomeryOne;
+        _previousExponent = 0UL;
+        _currentMontgomery = divisor.MontgomeryOne;
+        _hasState = false;
+    }
+
+    public bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL;
+
+    public bool HasState => _hasState;
+
+    public ulong PreviousExponent => _previousExponent;
+
+    public void Reset()
+    {
+        _previousExponent = 0UL;
+        _currentMontgomery = _montgomeryOne;
+        _hasState = false;
+    }
+
+    public ulong ComputeNext(ulong exponent)
+    {
+        if (!IsValidModulus)
+        {
+            _hasState = false;
+            return 0UL;
+        }
+
+        if (!_hasState || exponent <= _previousExponent)
+        {
+            _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor);
+            _previousExponent = exponent;
+            _hasState = true;
+            return ReduceCurrent();
+        }
+
+        ulong delta = exponent - _previousExponent;
+        ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor);
+        _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
+        _previousExponent = exponent;
+        return ReduceCurrent();
+    }
+
+    public bool ComputeNextIsUnity(ulong exponent)
+    {
+        if (!IsValidModulus)
+        {
+            _hasState = false;
+            return false;
+        }
+
+        if (!_hasState || exponent <= _previousExponent)
+        {
+            _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor);
+            _previousExponent = exponent;
+            _hasState = true;
+            return _currentMontgomery == _montgomeryOne;
+        }
+
+        ulong delta = exponent - _previousExponent;
+        ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor);
+        _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
+        _previousExponent = exponent;
+        return _currentMontgomery == _montgomeryOne;
+    }
+
+    public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity)
+    {
+        if (!IsValidModulus)
+        {
+            _hasState = false;
+            isUnity = false;
+            return false;
+        }
+
+        _currentMontgomery = montgomeryResult;
+        _previousExponent = exponent;
+        _hasState = true;
+        isUnity = _currentMontgomery == _montgomeryOne;
+        return true;
+    }
+
+    public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity)
+    {
+        if (!IsValidModulus)
+        {
+            _hasState = false;
+            isUnity = false;
+            return false;
+        }
+
+        if (!_hasState || exponent <= _previousExponent)
+        {
+            _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor);
+            _previousExponent = exponent;
+            _hasState = true;
+            isUnity = _currentMontgomery == _montgomeryOne;
+            return true;
+        }
+
+        _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime);
+        _previousExponent = exponent;
+        isUnity = _currentMontgomery == _montgomeryOne;
+        return true;
+    }
+
+    private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);
+}

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Buffers;
 using System.Collections.Concurrent;
 using System.Runtime.CompilerServices;
@@ -28,7 +29,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
 
         private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
-                _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputePrimeExponentKernel));
+                _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeMontgomeryExponentKernel));
 
 	public int GpuBatchSize
 	{
@@ -352,6 +353,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 private MemoryBuffer1D<ulong, Stride1D.Dense> _resultsBuffer;
                 private ulong[] _hostBuffer;
                 private int[] _positionBuffer;
+                private byte[] _modeBuffer;
                 private int _capacity;
                 private bool _disposed;
 
@@ -366,6 +368,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 			_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
                         _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
                         _positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
+                        _modeBuffer = ArrayPool<byte>.Shared.Rent(_capacity);
 		}
 
 		internal void Reset()
@@ -375,86 +378,128 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
                 public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
                 {
-			MersenneNumberDivisorByDivisorGpuTester owner = _owner;
+                        if (_disposed)
-			int gpuBatchSize = owner._gpuBatchSize,
+                        {
-				primesLength = primes.Length;
+                                throw new ObjectDisposedException(nameof(DivisorScanSession));
+                        }
 
-                        MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
+                        MersenneNumberDivisorByDivisorGpuTester owner = _owner;
-			ulong modulus = divisorData.Modulus;
+                        int gpuBatchSize = owner._gpuBatchSize;
-			if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+                        int primesLength = primes.Length;
+                        if (primesLength == 0)
                         {
                                 return;
                         }
 
+                        EnsureCapacity(gpuBatchSize);
+
+                        MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
+                        var exponentStepper = new ExponentRemainderStepper(divisorData);
+                        if (!exponentStepper.IsValidModulus)
+                        {
+                                hits.Clear();
+                                return;
+                        }
+
+                        bool useDivisorCycles = owner._useDivisorCycles;
+                        if (useDivisorCycles && divisorCycle == 0UL)
+                        {
+                                throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}.");
+                        }
+
+                        CycleRemainderStepper cycleStepper = useDivisorCycles ? new CycleRemainderStepper(divisorCycle) : default;
+
                         Accelerator accelerator = _accelerator;
                         Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel;
-			ArrayView1D<ulong, Stride1D.Dense> exponentView, exponentsView = _exponentsBuffer.View;
+                        ArrayView1D<ulong, Stride1D.Dense> exponentsView = _exponentsBuffer.View;
-			ArrayView1D<ulong, Stride1D.Dense> resultView, resultsView = _resultsBuffer.View;
+                        ArrayView1D<ulong, Stride1D.Dense> resultsView = _resultsBuffer.View;
 
-                        Span<ulong> exponentSlice, hostSpan = new(_hostBuffer);
+                        Span<ulong> hostSpan = _hostBuffer.AsSpan(0, _capacity);
-                        Span<int> positionSpan = new(_positionBuffer);
+                        Span<int> positionSpan = _positionBuffer.AsSpan(0, _capacity);
-                        Span<byte> hitsSlice;
+                        Span<byte> modeSpan = _modeBuffer.AsSpan(0, _capacity);
-                        ReadOnlySpan<ulong> primesSlice;
+
+                        int offset = 0;
 
-                        int batchSize, computeCount, i, offset = 0;
-                        bool useDivisorCycles = owner._useDivisorCycles && divisorCycle != 0UL;
-                        CycleRemainderStepper cycleStepper = default;
-                        if (useDivisorCycles)
-                        {
-                                cycleStepper = new CycleRemainderStepper(divisorCycle);
-                        }
                         while (offset < primesLength)
                         {
-                                batchSize = Math.Min(gpuBatchSize, primesLength - offset);
+                                int batchSize = Math.Min(gpuBatchSize, primesLength - offset);
-                                primesSlice = primes.Slice(offset, batchSize);
+                                ReadOnlySpan<ulong> primesSlice = primes.Slice(offset, batchSize);
-                                hitsSlice = hits.Slice(offset, batchSize);
+                                Span<byte> hitsSlice = hits.Slice(offset, batchSize);
+                                hitsSlice.Clear();
 
-                                computeCount = 0;
+                                int computeCount = 0;
+                                bool localHasState = exponentStepper.HasState;
+                                ulong localPreviousExponent = exponentStepper.PreviousExponent;
+
+                                for (int i = 0; i < batchSize; i++)
+                                {
+                                        ulong primeValue = primesSlice[i];
 
                                         if (useDivisorCycles)
                                         {
-                                        Span<ulong> remainderSpan = hostSpan.Slice(0, batchSize);
+                                                ulong remainder = cycleStepper.ComputeNext(primeValue);
-
-                                        for (i = 0; i < batchSize; i++)
-                                        {
-                                                ulong remainder = cycleStepper.ComputeNext(primesSlice[i]);
-                                                remainderSpan[i] = remainder;
-                                                hitsSlice[i] = 0;
-
                                                 if (remainder != 0UL)
                                                 {
                                                         continue;
                                                 }
-
-                                                positionSpan[computeCount] = i;
-                                                hostSpan[computeCount] = remainder;
-                                                computeCount++;
                                         }
+
+                                        ulong workValue;
+                                        byte mode;
+
+                                        if (!localHasState || primeValue <= localPreviousExponent)
+                                        {
+                                                mode = 0;
+                                                workValue = primeValue;
+                                                localHasState = true;
                                         }
                                         else
                                         {
-                                        for (i = 0; i < batchSize; i++)
+                                                mode = 1;
-                                        {
+                                                workValue = primeValue - localPreviousExponent;
-                                                positionSpan[computeCount] = i;
-                                                hostSpan[computeCount] = primesSlice[i];
-						computeCount++;
                                         }
+
+                                        localPreviousExponent = primeValue;
+
+                                        if (useDivisorCycles)
+                                        {
+                                                workValue %= divisorCycle;
+                                        }
+
+                                        positionSpan[computeCount] = i;
+                                        hostSpan[computeCount] = workValue;
+                                        modeSpan[computeCount] = mode;
+                                        computeCount++;
                                 }
 
                                 if (computeCount > 0)
                                 {
-                                        exponentSlice = hostSpan[..computeCount];
+                                        Span<ulong> exponentSlice = hostSpan[..computeCount];
-                                        exponentView = exponentsView.SubView(0, computeCount);
+                                        ArrayView1D<ulong, Stride1D.Dense> exponentView = exponentsView.SubView(0, computeCount);
-
                                         exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount);
-					resultView = resultsView.SubView(0, computeCount);
+
+                                        ArrayView1D<ulong, Stride1D.Dense> resultView = resultsView.SubView(0, computeCount);
                                         kernel(computeCount, divisorData, exponentView, resultView);
                                         accelerator.Synchronize();
                                         resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount);
 
-                                        for (i = 0; i < computeCount; i++)
+                                        for (int i = 0; i < computeCount; i++)
                                         {
-                                                hitsSlice[positionSpan[i]] = exponentSlice[i] == 1UL ? (byte)1 : (byte)0;
+                                                int position = positionSpan[i];
+                                                ulong primeValue = primesSlice[position];
+                                                ulong montgomeryResult = exponentSlice[i];
+                                                bool isUnity;
+
+                                                if (modeSpan[i] == 0)
+                                                {
+                                                        exponentStepper.TryInitializeFromMontgomeryResult(primeValue, montgomeryResult, out isUnity);
+                                                }
+                                                else
+                                                {
+                                                        exponentStepper.TryAdvanceWithMontgomeryDelta(primeValue, montgomeryResult, out isUnity);
+                                                }
+
+                                                hitsSlice[position] = isUnity ? (byte)1 : (byte)0;
                                         }
                                 }
 
@@ -462,6 +507,38 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         }
                 }
 
+                private void EnsureCapacity(int requiredCapacity)
+                {
+                        if (requiredCapacity <= _capacity)
+                        {
+                                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)
@@ -474,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 }
         }
 
-	private static void ComputePrimeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
+        private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
         {
                 ulong modulus = divisor.Modulus;
                 if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -483,7 +560,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         return;
                 }
 
-		results[index] = exponents[index].Pow2MontgomeryMod(divisor);
+                ulong exponent = exponents[index];
+
+                results[index] = exponent.Pow2MontgomeryModMontgomery(divisor);
         }
 
 	public ulong DivisorLimit

PerfectNumbers.Core/ULongExtensions.cs

@@ -246,6 +246,38 @@ public static class ULongExtensions
                 return result.MontgomeryMultiply(1UL, modulus, nPrime);
         }
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static ulong Pow2MontgomeryModMontgomery(this ulong exponent, in MontgomeryDivisorData divisor)
+        {
+                ulong modulus = divisor.Modulus;
+                if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+                {
+                        return 0UL;
+                }
+
+                ulong result = divisor.MontgomeryOne;
+                ulong baseVal = divisor.MontgomeryTwo;
+                ulong remainingExponent = exponent;
+
+                while (remainingExponent > 0UL)
+                {
+                        if ((remainingExponent & 1UL) != 0UL)
+                        {
+                                result = result.MontgomeryMultiply(baseVal, modulus, divisor.NPrime);
+                        }
+
+                        remainingExponent >>= 1;
+                        if (remainingExponent == 0UL)
+                        {
+                                break;
+                        }
+
+                        baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, divisor.NPrime);
+                }
+
+                return result;
+        }
+
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
         {