Merge pull request #90 from mlemanczyk/codex/optimize-memory-usage-in-evenperfectbitscanner

Optimize cycle remainder reuse in by-divisor testers

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

@@ -1,9 +1,9 @@
 using System;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
+using System.Numerics;
 using PerfectNumbers.Core;
 using PerfectNumbers.Core.Gpu;
-using System.Numerics;
 
 namespace EvenPerfectBitScanner.Benchmarks;
 
@@ -13,8 +13,8 @@ public class Pow2MontgomeryModBenchmarks
 {
     private const int SampleCount = 256;
 
-    private readonly MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[] _smallDivisors = new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[SampleCount];
+    private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
-    private readonly MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[] _largeDivisors = new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[SampleCount];
+    private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
     private readonly ulong[] _smallExponents = new ulong[SampleCount];
     private readonly ulong[] _largeExponents = new ulong[SampleCount];
     private readonly ulong[] _smallCycles = new ulong[SampleCount];
@@ -184,7 +184,7 @@ public class Pow2MontgomeryModBenchmarks
         return checksum;
     }
 
-    private void GetData(out ulong[] exponents, out MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[] divisors, out ulong[] cycles)
+    private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
     {
         if (Scale == InputScale.Small)
         {
@@ -217,14 +217,14 @@ public class Pow2MontgomeryModBenchmarks
         return (modulus, (ulong)bitLength);
     }
 
-    private static MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
+    private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
     {
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
         {
-            return new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
+            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
         }
 
-        return new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData(
+        return new MontgomeryDivisorData(
             modulus,
             ComputeMontgomeryNPrime(modulus),
             ComputeMontgomeryResidue(1UL, modulus),
@@ -245,7 +245,7 @@ public class Pow2MontgomeryModBenchmarks
         return unchecked(0UL - inv);
     }
 
-    private static ulong Pow2MontgomeryModLeftToRight(ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor)
+    private static ulong Pow2MontgomeryModLeftToRight(ulong exponent, in MontgomeryDivisorData divisor)
     {
         ulong modulus = divisor.Modulus;
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -277,7 +277,7 @@ public class Pow2MontgomeryModBenchmarks
         return result.MontgomeryMultiply(1UL, modulus, nPrime);
     }
 
-    private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor)
+    private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor)
     {
         ulong modulus = divisor.Modulus;
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -326,7 +326,7 @@ public class Pow2MontgomeryModBenchmarks
         return result.MontgomeryMultiply(1UL, modulus, nPrime);
     }
 
-    private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor)
+    private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor)
     {
         ulong modulus = divisor.Modulus;
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)

EvenPerfectBitScanner/Program.cs

@@ -899,9 +899,16 @@ internal static class Program
 							divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining);
 							if (divisor != 0UL)
 							{
-								divisorCycle = useDivisorCycles
+							DivisorCycleCache.CycleBlock? cycleBlock = null;
-										? DivisorCycleCache.Shared.Acquire(divisor).GetCycle(divisor)
+							if (useDivisorCycles)
-										: 0UL;
+							{
+							cycleBlock = DivisorCycleCache.Shared.Acquire(divisor);
+							divisorCycle = cycleBlock.GetCycle(divisor);
+							}
+							else
+							{
+							divisorCycle = 0UL;
+							}
 
 							activeCount = BuildPrimeBuffer(divisor, primeValues, allowedMaxValues, stateFlags, primeBuffer, indexBuffer, completionsBuffer, ref completionsCount, ref remainingStates, activeStateMask, ref activeStartIndex);
 
@@ -912,6 +919,7 @@ internal static class Program
 
 							if (activeCount == 0)
 							{
+							cycleBlock?.Dispose();
 							continue;
 							}
 
@@ -945,6 +953,8 @@ internal static class Program
 							FlushPendingResults(compositesBuffer, ref compositesCount);
 							}
 
+							cycleBlock?.Dispose();
+
 								continue;
 							}
 

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -1,6 +1,7 @@
 using FluentAssertions;
 using PerfectNumbers.Core;
 using PerfectNumbers.Core.Gpu;
+using PerfectNumbers.Core.Cpu;
 using System.Reflection;
 using Xunit;
 
@@ -37,7 +38,7 @@ public class MersenneNumberDivisorGpuTesterTests
         var tester = new MersenneNumberDivisorGpuTester();
         typeof(MersenneNumberDivisorGpuTester)
             .GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)!
-            .SetValue(null, Array.Empty<(ulong, ulong)>());
+            .SetValue(null, Array.Empty<(ulong, uint)>());
 
         tester.IsPrime(11UL, UInt128.Zero, 0UL, out bool divisorsExhausted).Should().BeTrue();
         divisorsExhausted.Should().BeFalse();
@@ -102,6 +103,22 @@ public class MersenneNumberDivisorGpuTesterTests
         divisorsExhausted.Should().BeTrue();
     }
 
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_gpu_tester_uses_cycle_remainders_per_divisor()
+    {
+        var tester = new MersenneNumberDivisorByDivisorGpuTester
+        {
+            UseDivisorCycles = true,
+            GpuBatchSize = 8,
+        };
+
+        tester.ConfigureFromMaxPrime(13UL);
+
+        tester.IsPrime(11UL, out bool divisorsExhausted).Should().BeFalse();
+        divisorsExhausted.Should().BeTrue();
+    }
+
     [Fact]
     [Trait("Category", "Fast")]
     public void ByDivisor_session_checks_divisors_across_primes()
@@ -122,5 +139,49 @@ public class MersenneNumberDivisorGpuTesterTests
         session.CheckDivisor(31UL, cycle31, primes, hits);
         hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
     }
+
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_gpu_session_reuses_cycle_remainders_across_batches()
+    {
+        var tester = new MersenneNumberDivisorByDivisorGpuTester
+        {
+            UseDivisorCycles = true,
+            GpuBatchSize = 2,
+        };
+
+        tester.ConfigureFromMaxPrime(17UL);
+
+        using var session = tester.CreateDivisorSession();
+        ulong[] primes = { 5UL, 7UL, 11UL, 13UL, 17UL };
+        byte[] hits = new byte[primes.Length];
+
+        ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL);
+        session.CheckDivisor(23UL, cycle23, primes, hits);
+
+        hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0, 0 });
+    }
+
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_cpu_session_reuses_cycle_remainders_across_primes()
+    {
+        var tester = new MersenneNumberDivisorByDivisorCpuTester
+        {
+            UseDivisorCycles = true,
+            BatchSize = 2,
+        };
+
+        tester.ConfigureFromMaxPrime(17UL);
+
+        using var session = tester.CreateDivisorSession();
+        ulong[] primes = { 5UL, 7UL, 11UL, 13UL, 17UL };
+        byte[] hits = new byte[primes.Length];
+
+        ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL);
+        session.CheckDivisor(23UL, cycle23, primes, hits);
+
+        hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0, 0 });
+    }
 }
 

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -1,9 +1,7 @@
 using System;
 using System.Collections.Concurrent;
-using System.Runtime.CompilerServices;
+using System.Numerics;
 using PerfectNumbers.Core;
-using PerfectNumbers.Core.Gpu;
-using MontgomeryDivisorData = PerfectNumbers.Core.Gpu.MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData;
 
 namespace PerfectNumbers.Core.Cpu;
 
@@ -167,12 +165,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             processedCount++;
             lastProcessed = divisor;
 
-            MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
+            MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
             byte hit;
+            ulong divisorCycle = 0UL;
+            DivisorCycleCache.CycleBlock? cycleBlock = null;
             if (useCycles)
             {
-                DivisorCycleCache.CycleBlock cycleBlock = DivisorCycleCache.Shared.Acquire(divisor);
+                cycleBlock = DivisorCycleCache.Shared.Acquire(divisor);
-                ulong divisorCycle = cycleBlock.GetCycle(divisor);
+                divisorCycle = cycleBlock.GetCycle(divisor);
                 hit = CheckDivisor(prime, divisorCycle != 0UL, divisorCycle, divisorData);
             }
             else
@@ -180,6 +180,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 hit = CheckDivisor(prime, false, 0UL, divisorData);
             }
 
+            cycleBlock?.Dispose();
+
             if (hit != 0)
             {
                 processedAll = true;
@@ -270,35 +272,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
     }
 
-    private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
-    {
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
-        }
-
-        return new MontgomeryDivisorData(
-                modulus,
-                ComputeMontgomeryNPrime(modulus),
-                ComputeMontgomeryResidue(1UL, modulus),
-                ComputeMontgomeryResidue(2UL, modulus));
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
-
-    private static ulong ComputeMontgomeryNPrime(ulong modulus)
-    {
-        ulong inv = modulus;
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        return unchecked(0UL - inv);
-    }
-
     private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
     {
         private readonly MersenneNumberDivisorByDivisorCpuTester _owner;
@@ -335,26 +308,31 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             }
 
             bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL;
-            MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
+            MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
-
-            for (int i = 0; i < length; i++)
-            {
-                ulong prime = primes[i];
 
             if (cycleEnabled)
             {
-                    ulong remainder = prime % divisorCycle;
+                var stepper = new CycleRemainderStepper(divisorCycle);
+
+                for (int i = 0; i < length; i++)
+                {
+                    ulong remainder = stepper.ComputeNext(primes[i]);
                     if (remainder != 0UL)
                     {
                         hits[i] = 0;
                         continue;
                     }
+
+                    ulong residue = remainder.Pow2MontgomeryModFromCycleRemainder(divisorData);
+                    hits[i] = residue == 1UL ? (byte)1 : (byte)0;
                 }
 
-                ulong residue = cycleEnabled
+                return;
-                    ? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData)
+            }
-                    : prime.Pow2MontgomeryMod(divisorData);
 
+            for (int i = 0; i < length; i++)
+            {
+                ulong residue = primes[i].Pow2MontgomeryMod(divisorData);
                 hits[i] = residue == 1UL ? (byte)1 : (byte)0;
             }
         }

PerfectNumbers.Core/CycleRemainderStepper.cs

@@ -0,0 +1,84 @@
+using System.Numerics;
+
+namespace PerfectNumbers.Core;
+
+internal struct CycleRemainderStepper
+{
+    private readonly ulong _cycleLength;
+    private ulong _previousPrime;
+    private ulong _currentRemainder;
+    private bool _hasState;
+
+    public CycleRemainderStepper(ulong cycleLength)
+    {
+        _cycleLength = cycleLength;
+        _previousPrime = 0UL;
+        _currentRemainder = 0UL;
+        _hasState = false;
+    }
+
+    public bool HasCycle => _cycleLength != 0UL;
+
+    public void Reset()
+    {
+        _previousPrime = 0UL;
+        _currentRemainder = 0UL;
+        _hasState = false;
+    }
+
+    public void InitializeState(ulong previousPrime, ulong previousRemainder, bool hasState)
+    {
+        _previousPrime = previousPrime;
+        _currentRemainder = previousRemainder;
+        _hasState = hasState;
+    }
+
+    public ulong ComputeNext(ulong prime)
+    {
+        if (_cycleLength == 0UL)
+        {
+            _previousPrime = prime;
+            _currentRemainder = 0UL;
+            _hasState = true;
+            return 0UL;
+        }
+
+        if (!_hasState || prime <= _previousPrime)
+        {
+            _currentRemainder = prime % _cycleLength;
+            _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;
+                    }
+                }
+            }
+        }
+
+        _previousPrime = prime;
+        return _currentRemainder;
+    }
+
+    public void CaptureState(out ulong previousPrime, out ulong previousRemainder, out bool hasState)
+    {
+        previousPrime = _previousPrime;
+        previousRemainder = _currentRemainder;
+        hasState = _hasState;
+    }
+}

PerfectNumbers.Core/DivisorCycleCache.cs

@@ -3,6 +3,7 @@ using System.Buffers;
 using System.Collections.Concurrent;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using System.Threading;
 using System.Threading.Tasks;
 using ILGPU;
 using ILGPU.Runtime;
@@ -12,15 +13,51 @@ namespace PerfectNumbers.Core;
 
 public sealed class DivisorCycleCache
 {
-    public sealed class CycleBlock(int index, ulong start, ulong[] cycles)
+    public sealed class CycleBlock : IDisposable
     {
-		internal readonly int Index = index;
+        private readonly DivisorCycleCache _owner;
-		internal readonly ulong Start = start;
+        private int _referenceCount;
-		internal readonly ulong End = start + (ulong)cycles.Length - 1UL;
+
-		internal readonly ulong[] Cycles = cycles;
+        internal CycleBlock(DivisorCycleCache owner, int index, ulong start, in ulong[] cycles)
+        {
+            _owner = owner;
+            Index = index;
+            Start = start;
+            End = start + (ulong)cycles.Length - 1UL;
+            Cycles = cycles;
+        }
+
+        internal int Index { get; }
+
+        internal ulong Start { get; }
+
+        internal ulong End { get; }
+
+        internal ulong[] Cycles { get; }
+
+        internal int ReferenceCount => Volatile.Read(ref _referenceCount);
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public ulong GetCycle(ulong divisor) => Cycles[(int)(divisor - Start)];
+
+        internal void Retain()
+        {
+            Interlocked.Increment(ref _referenceCount);
+        }
+
+        public void Dispose()
+        {
+            int remaining = Interlocked.Decrement(ref _referenceCount);
+            if (remaining < 0)
+            {
+                throw new InvalidOperationException("Cycle block released more times than it was acquired.");
+            }
+
+            if (remaining == 0)
+            {
+                _owner.OnCycleBlockReleased(this);
+            }
+        }
     }
 
     private const int CycleGenerationBatchSize = 262_144;
@@ -63,24 +100,37 @@ public sealed class DivisorCycleCache
 
     private void Initialize(ulong[] snapshot)
     {
+        CycleBlock? activeToDispose = null;
+        CycleBlock? prefetchedToDispose = null;
+
         lock (_sync)
         {
+            if (_activeBlock is not null && _activeBlock.Index != 0)
+            {
+                activeToDispose = _activeBlock;
+            }
+
+            if (_prefetchedBlock is not null && _prefetchedBlock.Index != 0)
+            {
+                prefetchedToDispose = _prefetchedBlock;
+            }
+
             _pending.Clear();
-            _baseBlock = new CycleBlock(index: 0, start: 0UL, snapshot);
+            _baseBlock = new CycleBlock(this, index: 0, start: 0UL, snapshot);
             _activeBlock = null;
             _prefetchedBlock = null;
             StartPrefetchLocked(1);
         }
+
+        activeToDispose?.Dispose();
+        prefetchedToDispose?.Dispose();
     }
 
     public CycleBlock Acquire(ulong divisor)
     {
-        if (divisor <= _baseBlock.End)
+        CycleBlock block = divisor <= _baseBlock.End ? _baseBlock : AcquireDynamicBlock(divisor);
-        {
+        block.Retain();
-            return _baseBlock;
+        return block;
-        }
-
-        return AcquireDynamicBlock(divisor);
     }
 
     private CycleBlock AcquireDynamicBlock(ulong divisor)
@@ -149,21 +199,42 @@ public sealed class DivisorCycleCache
             return;
         }
 
+        CycleBlock? previousActive = _activeBlock;
         _activeBlock = prefetched;
         _prefetchedBlock = null;
+        prefetched.Retain();
+        prefetched.Dispose();
+
+        if (previousActive is not null && previousActive.Index != 0 && previousActive.Index != prefetched.Index)
+        {
+            previousActive.Dispose();
+        }
+
         StartPrefetchLocked(prefetched.Index + 1);
     }
 
     private void SetActiveBlockLocked(CycleBlock block)
     {
+        CycleBlock? previousActive = _activeBlock;
         _activeBlock = block;
+        block.Retain();
 
         CycleBlock? prefetched = _prefetchedBlock;
         if (prefetched is not null && prefetched.Index != block.Index + 1)
         {
+            if (prefetched.Index != 0)
+            {
+                prefetched.Dispose();
+            }
+
             _prefetchedBlock = null;
         }
 
+        if (previousActive is not null && previousActive.Index != 0 && previousActive.Index != block.Index)
+        {
+            previousActive.Dispose();
+        }
+
         StartPrefetchLocked(block.Index + 1);
     }
 
@@ -216,12 +287,59 @@ public sealed class DivisorCycleCache
         CycleBlock? active = _activeBlock;
         if (active is not null && blockIndex != active.Index + 1)
         {
+            MontgomeryDivisorDataCache.ReleaseBlock(block);
             return;
         }
 
-        if (_prefetchedBlock is null || _prefetchedBlock.Index < blockIndex)
+        CycleBlock? previousPrefetched = _prefetchedBlock;
+        if (ReferenceEquals(previousPrefetched, block))
         {
+            return;
+        }
+
+        if (_prefetchedBlock is null || _prefetchedBlock.Index <= blockIndex)
+        {
+            if (previousPrefetched is not null && previousPrefetched.Index != 0)
+            {
+                previousPrefetched.Dispose();
+            }
+
             _prefetchedBlock = block;
+            block.Retain();
+        }
+        else
+        {
+            MontgomeryDivisorDataCache.ReleaseBlock(block);
+        }
+    }
+
+    private void OnCycleBlockReleased(CycleBlock block)
+    {
+        if (block.Index == 0)
+        {
+            return;
+        }
+
+        for (int attempt = 0; attempt < 8; attempt++)
+        {
+            if (block.ReferenceCount != 0)
+            {
+                return;
+            }
+
+            Thread.SpinWait(1 << attempt);
+        }
+
+        bool stillTracked;
+
+        lock (_sync)
+        {
+            stillTracked = ReferenceEquals(_activeBlock, block) || ReferenceEquals(_prefetchedBlock, block);
+        }
+
+        if (!stillTracked && block.ReferenceCount == 0)
+        {
+            MontgomeryDivisorDataCache.ReleaseBlock(block);
         }
     }
 
@@ -245,7 +363,7 @@ public sealed class DivisorCycleCache
             ComputeCyclesCpu(start, cycles);
         }
 
-        return new CycleBlock(blockIndex, start, cycles);
+        return new CycleBlock(this, blockIndex, start, cycles);
     }
 
     private int GetBlockIndex(ulong divisor)

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -2,8 +2,12 @@ using System.Buffers;
 using System.Collections.Concurrent;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using System.Numerics;
 using ILGPU;
 using ILGPU.Runtime;
+using MontgomeryDivisorData = PerfectNumbers.Core.MontgomeryDivisorData;
+using MontgomeryDivisorDataCache = PerfectNumbers.Core.MontgomeryDivisorDataCache;
+using CycleRemainderStepper = PerfectNumbers.Core.CycleRemainderStepper;
 
 namespace PerfectNumbers.Core.Gpu;
 
@@ -15,21 +19,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         private bool _isConfigured;
         private bool _useDivisorCycles;
 
-	public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo)
+        private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>>> _kernelCache = new();
-	{
-		public readonly ulong Modulus = modulus;
-		public readonly ulong NPrime = nPrime;
-		public readonly ulong MontgomeryOne = montgomeryOne;
-		public readonly ulong MontgomeryTwo = montgomeryTwo;
-	}
-
-	private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, byte, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>> _kernelCache = new();
 	private readonly ConcurrentDictionary<Accelerator, Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>> _kernelExponentCache = new();
 	private readonly ConcurrentDictionary<Accelerator, ConcurrentBag<BatchResources>> _resourcePools = new(AcceleratorReferenceComparer.Instance);
 	private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
 
-	private Action<Index1D, ulong, byte, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
+        private Action<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
-																	_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, byte, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>(CheckKernel));
+                _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));
 
@@ -103,8 +100,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                                                         accelerator,
                                                         kernel,
                                                         resources.DivisorsBuffer,
+                                                        resources.ExponentBuffer,
                                                         resources.HitsBuffer,
                                                         resources.Divisors,
+                                                        resources.Exponents,
                                                         resources.Hits,
                                                         resources.DivisorData,
                                                         out lastProcessed,
@@ -133,10 +132,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         ulong allowedMax,
                         bool useCycles,
                         Accelerator accelerator,
-			Action<Index1D, ulong, byte, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> kernel,
+                        Action<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>> kernel,
                         MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> divisorsBuffer,
+                        MemoryBuffer1D<ulong, Stride1D.Dense> exponentBuffer,
                         MemoryBuffer1D<byte, Stride1D.Dense> hitsBuffer,
                         in ulong[] divisors,
+                        ulong[] exponents,
                         byte[] hits,
                         MontgomeryDivisorData[] divisorData,
 			out ulong lastProcessed,
@@ -156,14 +157,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
                 bool composite = false;
                 bool processedAll = false;
-		ulong currentDivisor = 3UL, cycle = 0UL, next, nextDivisor;
+                ulong currentDivisor = 3UL, next, nextDivisor;
                 int batchSize, i;
                 bool reachedEndInBatch;
 
                 Span<MontgomeryDivisorData> divisorDataSpan;
                 Span<ulong> divisorSpan;
+                Span<ulong> exponentSpan;
                 Span<byte> hitsSpan;
                 ArrayView1D<MontgomeryDivisorData, Stride1D.Dense> divisorView;
+                ArrayView1D<ulong, Stride1D.Dense> exponentView;
                 ArrayView1D<byte, Stride1D.Dense> hitsView;
                 DivisorCycleCache divisorCyclesCache = DivisorCycleCache.Shared;
                 DivisorCycleCache.CycleBlock? cycleLease = useCycles ? divisorCyclesCache.Acquire(currentDivisor) : null;
@@ -180,24 +183,37 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 				processedCount++;
 
                                 bool includeDivisor = true;
+                                ulong exponent = prime;
+
                                 if (useCycles)
                                 {
                                         if (currentDivisor > cycleLeaseEnd)
                                         {
+                                                cycleLease?.Dispose();
                                                 cycleLease = divisorCyclesCache.Acquire(currentDivisor);
                                                 cycleLeaseEnd = cycleLease.End;
                                         }
 
-					cycle = prime % cycleLease!.GetCycle(currentDivisor);
+                                        ulong cycleLength = cycleLease!.GetCycle(currentDivisor);
-					if (cycle != 0UL)
+                                        if (cycleLength > 0UL)
+                                        {
+                                                ulong remainder = prime % cycleLength;
+                                                if (remainder != 0UL)
                                                 {
                                                         includeDivisor = false;
                                                 }
+                                                else
+                                                {
+                                                        exponent = remainder;
+                                                }
+                                        }
                                 }
 
                                 if (includeDivisor)
                                 {
-					divisors[batchSize++] = currentDivisor;
+                                        divisors[batchSize] = currentDivisor;
+                                        exponents[batchSize] = exponent;
+                                        batchSize++;
                                 }
 
 				lastProcessed = currentDivisor;
@@ -232,19 +248,21 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
                         divisorDataSpan = divisorData.AsSpan(0, batchSize);
                         divisorView = divisorsBuffer.View.SubView(0, batchSize);
+                        exponentView = exponentBuffer.View.SubView(0, batchSize);
                         hitsView = hitsBuffer.View.SubView(0, batchSize);
 
                         divisorSpan = divisors.AsSpan(0, batchSize);
+                        exponentSpan = exponents.AsSpan(0, batchSize);
                         hitsSpan = hits.AsSpan(0, batchSize);
 
                         for (i = 0; i < batchSize; i++)
                         {
-				divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]);
+                                divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]);
                         }
 
                         divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize);
-			byte useCycleChecks = useCycles ? (byte)1 : (byte)0;
+                        exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), batchSize);
-			kernel(batchSize, prime, useCycleChecks, cycle, divisorView, hitsView);
+                        kernel(batchSize, divisorView, exponentView, hitsView);
 			accelerator.Synchronize();
 			hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize);
 
@@ -274,6 +292,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 			}
 		}
 
+                cycleLease?.Dispose();
                 coveredRange = composite || processedAll || currentDivisor > allowedMax;
 
                 return composite;
@@ -309,7 +328,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 		return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
 	}
 
-	private static void CheckKernel(Index1D index, ulong prime, byte useCycleChecks, ulong cycle, ArrayView<MontgomeryDivisorData> divisors, ArrayView<byte> hits)
+        private static void CheckKernel(Index1D index, ArrayView<MontgomeryDivisorData> divisors, ArrayView<ulong> exponents, ArrayView<byte> hits)
         {
                 MontgomeryDivisorData divisor = divisors[index];
                 ulong modulus = divisor.Modulus;
@@ -319,12 +338,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         return;
                 }
 
-		if (useCycleChecks != 0)
+                ulong exponent = exponents[index];
-		{
+                hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0;
-			prime = cycle;
-		}
-
-		hits[index] = prime.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0;
         }
 
 	public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -364,7 +379,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 			int gpuBatchSize = owner._gpuBatchSize,
 				primesLength = primes.Length;
 
-			MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
+                        MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
 			ulong modulus = divisorData.Modulus;
 			if (modulus <= 1UL || (modulus & 1UL) == 0UL)
 			{
@@ -382,8 +397,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         ReadOnlySpan<ulong> primesSlice;
 
                         int batchSize, computeCount, i, offset = 0;
-			bool useDivisorCycles = owner._useDivisorCycles;
+                        bool useDivisorCycles = owner._useDivisorCycles && divisorCycle != 0UL;
-			ulong residue;
+                        CycleRemainderStepper cycleStepper = default;
+                        if (useDivisorCycles)
+                        {
+                                cycleStepper = new CycleRemainderStepper(divisorCycle);
+                        }
                         while (offset < primesLength)
                         {
                                 batchSize = Math.Min(gpuBatchSize, primesLength - offset);
@@ -394,17 +413,21 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
                                 if (useDivisorCycles)
                                 {
+                                        Span<ulong> remainderSpan = hostSpan.Slice(0, batchSize);
+
                                         for (i = 0; i < batchSize; i++)
                                         {
-						residue = primesSlice[i] % divisorCycle;
+                                                ulong remainder = cycleStepper.ComputeNext(primesSlice[i]);
-						if (residue != 0UL)
+                                                remainderSpan[i] = remainder;
+                                                hitsSlice[i] = 0;
+
+                                                if (remainder != 0UL)
                                                 {
-							// We don't need to clear hitsSlice[i] because the parent function calls .Clear() on the buffer
                                                         continue;
                                                 }
 
                                                 positionSpan[computeCount] = i;
-						hostSpan[computeCount] = residue;
+                                                hostSpan[computeCount] = remainder;
                                                 computeCount++;
                                         }
                                 }
@@ -463,35 +486,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 		results[index] = exponents[index].Pow2MontgomeryMod(divisor);
 	}
 
-	private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
-	{
-		if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-		{
-			return new(modulus, 0UL, 0UL, 0UL);
-		}
-
-		return new(
-			modulus,
-			ComputeMontgomeryNPrime(modulus),
-			ComputeMontgomeryResidue(1UL, modulus),
-			ComputeMontgomeryResidue(2UL, modulus)
-		);
-	}
-
-	private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
-
-	private static ulong ComputeMontgomeryNPrime(ulong modulus)
-	{
-		ulong inv = modulus;
-		inv *= unchecked(2UL - modulus * inv);
-		inv *= unchecked(2UL - modulus * inv);
-		inv *= unchecked(2UL - modulus * inv);
-		inv *= unchecked(2UL - modulus * inv);
-		inv *= unchecked(2UL - modulus * inv);
-		inv *= unchecked(2UL - modulus * inv);
-		return unchecked(0UL - inv);
-	}
-
 	public ulong DivisorLimit
 	{
 		get
@@ -526,8 +520,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 internal BatchResources(Accelerator accelerator, int capacity)
                 {
                         DivisorsBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(capacity);
+                        ExponentBuffer = accelerator.Allocate1D<ulong>(capacity);
                         HitsBuffer = accelerator.Allocate1D<byte>(capacity);
                         Divisors = ArrayPool<ulong>.Shared.Rent(capacity);
+                        Exponents = ArrayPool<ulong>.Shared.Rent(capacity);
                         Hits = ArrayPool<byte>.Shared.Rent(capacity);
                         DivisorData = ArrayPool<MontgomeryDivisorData>.Shared.Rent(capacity);
                         Capacity = capacity;
@@ -535,10 +531,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
                 internal readonly MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> DivisorsBuffer;
 
+                internal readonly MemoryBuffer1D<ulong, Stride1D.Dense> ExponentBuffer;
+
                 internal readonly MemoryBuffer1D<byte, Stride1D.Dense> HitsBuffer;
 
                 internal readonly ulong[] Divisors;
 
+                internal readonly ulong[] Exponents;
+
                 internal readonly byte[] Hits;
 
                 internal readonly MontgomeryDivisorData[] DivisorData;
@@ -548,8 +548,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 public void Dispose()
                 {
                         DivisorsBuffer.Dispose();
+                        ExponentBuffer.Dispose();
                         HitsBuffer.Dispose();
                         ArrayPool<ulong>.Shared.Return(Divisors, clearArray: false);
+                        ArrayPool<ulong>.Shared.Return(Exponents, clearArray: false);
                         ArrayPool<byte>.Shared.Return(Hits, clearArray: false);
                         ArrayPool<MontgomeryDivisorData>.Shared.Return(DivisorData, clearArray: false);
                 }

PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs

@@ -16,7 +16,7 @@ public sealed class MersenneNumberDivisorGpuTester
         public static void BuildDivisorCandidates()
         {
                 ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
-                (ulong divisor, ulong cycle)[] list = new (ulong divisor, ulong cycle)[snapshot.Length / 2];
+                (ulong divisor, uint cycle)[] list = new (ulong divisor, uint cycle)[snapshot.Length / 2];
                 ulong cycle;
                 int count = 0, i, snapshotLength = snapshot.Length;
                 for (i = 3; i < snapshotLength; i += 2)
@@ -27,7 +27,7 @@ public sealed class MersenneNumberDivisorGpuTester
                                 continue;
                         }
 
-			list[count++] = ((ulong)i, cycle);
+                        list[count++] = ((ulong)i, (uint)cycle);
                 }
 
                 _divisorCandidates = count == 0 ? [] : list[..count];
@@ -52,7 +52,7 @@ public sealed class MersenneNumberDivisorGpuTester
 		return divisible;
 	}
 
-        private static (ulong divisor, ulong cycle)[]? _divisorCandidates = Array.Empty<(ulong divisor, ulong cycle)>();
+        private static (ulong divisor, uint cycle)[]? _divisorCandidates = Array.Empty<(ulong divisor, uint cycle)>();
 
 	public bool IsPrime(ulong p, UInt128 d, ulong divisorCyclesSearchLimit, out bool divisorsExhausted)
 	{
@@ -71,7 +71,7 @@ public sealed class MersenneNumberDivisorGpuTester
                 if (_divisorCandidates is { Length: > 0 } candidates)
                 {
                         int len = candidates.Length;
-                        ulong cycle;
+                        uint cycle;
                         for (int k = 0; k < len; k++)
                         {
                                 (ulong dSmall, cycle) = candidates[k];

PerfectNumbers.Core/MontgomeryDivisorData.cs

@@ -0,0 +1,148 @@
+using System.Collections.Concurrent;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo)
+{
+    public ulong Modulus { get; } = modulus;
+
+    public ulong NPrime { get; } = nPrime;
+
+    public ulong MontgomeryOne { get; } = montgomeryOne;
+
+    public ulong MontgomeryTwo { get; } = montgomeryTwo;
+}
+
+internal static class MontgomeryDivisorDataCache
+{
+    private const ulong BaseBlockEnd = PerfectNumberConstants.MaxQForDivisorCycles;
+
+    private static readonly ConcurrentDictionary<ulong, CacheEntry> Cache = new();
+    private static readonly ConcurrentDictionary<int, ConcurrentDictionary<ulong, byte>> BlockMembership = new();
+
+    public static MontgomeryDivisorData Get(ulong modulus)
+    {
+        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        {
+            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
+        }
+
+        int blockIndex = ComputeBlockIndex(modulus);
+        return GetManaged(modulus, blockIndex);
+    }
+
+    internal static void ReleaseBlock(DivisorCycleCache.CycleBlock block)
+    {
+        if (block.Index <= 0)
+        {
+            return;
+        }
+
+        if (BlockMembership.TryRemove(block.Index, out ConcurrentDictionary<ulong, byte>? members))
+        {
+            foreach (ulong modulus in members.Keys)
+            {
+                Cache.TryRemove(modulus, out _);
+            }
+        }
+    }
+
+    private static MontgomeryDivisorData GetManaged(ulong modulus, int blockIndex)
+    {
+        while (true)
+        {
+            if (Cache.TryGetValue(modulus, out CacheEntry entry))
+            {
+                if (entry.BlockIndex != blockIndex)
+                {
+                    if (blockIndex > 0 && entry.BlockIndex <= 0)
+                    {
+                        CacheEntry updated = new(entry.Data, blockIndex);
+                        if (Cache.TryUpdate(modulus, updated, entry))
+                        {
+                            RegisterBlockMembership(modulus, blockIndex);
+                            return entry.Data;
+                        }
+
+                        continue;
+                    }
+
+                    if (blockIndex > 0)
+                    {
+                        RegisterBlockMembership(modulus, blockIndex);
+                    }
+
+                    return entry.Data;
+                }
+                return entry.Data;
+            }
+
+            MontgomeryDivisorData created = Create(modulus);
+            if (Cache.TryAdd(modulus, new CacheEntry(created, blockIndex)))
+            {
+                RegisterBlockMembership(modulus, blockIndex);
+                return created;
+            }
+        }
+    }
+
+    private static void RegisterBlockMembership(ulong modulus, int blockIndex)
+    {
+        if (blockIndex <= 0)
+        {
+            return;
+        }
+
+        ConcurrentDictionary<ulong, byte> set = BlockMembership.GetOrAdd(blockIndex, static _ => new ConcurrentDictionary<ulong, byte>());
+        set.TryAdd(modulus, 0);
+    }
+
+    private static int ComputeBlockIndex(ulong modulus)
+    {
+        if (modulus <= BaseBlockEnd)
+        {
+            return 0;
+        }
+
+        return (int)((modulus - (BaseBlockEnd + 1UL)) / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1;
+    }
+
+    private static MontgomeryDivisorData Create(ulong modulus)
+    {
+        return new MontgomeryDivisorData(
+            modulus,
+            ComputeMontgomeryNPrime(modulus),
+            ComputeMontgomeryResidue(1UL, modulus),
+            ComputeMontgomeryResidue(2UL, modulus));
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
+
+    private static ulong ComputeMontgomeryNPrime(ulong modulus)
+    {
+        ulong inv = modulus;
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        inv *= unchecked(2UL - modulus * inv);
+        return unchecked(0UL - inv);
+    }
+
+    private readonly struct CacheEntry
+    {
+        internal CacheEntry(MontgomeryDivisorData data, int blockIndex)
+        {
+            Data = data;
+            BlockIndex = blockIndex;
+        }
+
+        internal MontgomeryDivisorData Data { get; }
+
+        internal int BlockIndex { get; }
+    }
+}

PerfectNumbers.Core/ULongExtensions.cs

@@ -214,7 +214,7 @@ public static class ULongExtensions
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static ulong Pow2MontgomeryMod(this ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor)
+        public static ulong Pow2MontgomeryMod(this ulong exponent, in MontgomeryDivisorData divisor)
         {
                 ulong modulus = divisor.Modulus;
                 if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -247,7 +247,7 @@ public static class ULongExtensions
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor)
+        public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
         {
                 if (cycleLength == 0UL)
                 {
@@ -255,12 +255,18 @@ public static class ULongExtensions
                 }
 
                 ulong rotationCount = exponent % cycleLength;
-                if (rotationCount == 0UL)
+                return rotationCount.Pow2MontgomeryModFromCycleRemainder(divisor);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+        {
+                if (reducedExponent == 0UL)
                 {
                         return 1UL;
                 }
 
-                return rotationCount.Pow2MontgomeryMod(divisor);
+                return reducedExponent.Pow2MontgomeryMod(divisor);
         }
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]