Document by-divisor option in CLI help

EvenPerfectBitScanner/Program.cs

@@ -30,8 +30,10 @@ internal static class Program
 	private static bool _useGcdFilter;
         private static bool _useDivisor;
         private static bool _useResidueMode;
+        private static bool _useByDivisorMode;
         private static UInt128 _divisor;
         private static MersenneNumberDivisorGpuTester? _divisorTester;
+        private static MersenneNumberDivisorByDivisorGpuTester? _byDivisorTester;
 	private static ulong? _orderWarmupLimitOverride;
 	private static unsafe delegate*<ulong, ref ulong, ulong> _transformP;
 	private static long _state;
@@ -68,6 +70,7 @@ internal static class Program
 		bool useLucas = false;
                 bool useResidue = false;    // M_p test via residue divisors (replaces LL/incremental)
                 bool useDivisor = false;     // M_p divisibility by specific divisor
+                bool useByDivisor = false;   // Iterative divisor scan across primes
 		UInt128 divisor = UInt128.Zero;
 		// Device routing
 		bool useGpuCycles = true;
@@ -131,6 +134,13 @@ internal static class Program
                                         useLucas = false;
                                         useResidue = false;
                                 }
+                                else if (value.Equals("bydivisor", StringComparison.OrdinalIgnoreCase))
+                                {
+                                        useByDivisor = true;
+                                        useLucas = false;
+                                        useResidue = false;
+                                        useDivisor = false;
+                                }
 				else
 				{
 					kernelType = GpuKernelType.Incremental;
@@ -419,8 +429,15 @@ internal static class Program
 
 		Console.WriteLine("Divisor cycles are ready");
 
+                if (useByDivisor)
+                {
+                        threadCount = 1;
+                        blockSize = 1;
+                }
+
                 _useDivisor = useDivisor;
                 _useResidueMode = useResidue;
+                _useByDivisorMode = useByDivisor;
                 _divisor = divisor;
 		if (useBitTransform)
 		{
@@ -438,7 +455,7 @@ internal static class Program
 		PrimeTesters = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
 		// Note: --primes-device controls default device for library kernels; p primality remains CPU here.
 		// Initialize per-thread p residue tracker (Identity model) at currentP
-		if (!useDivisor)
+                if (!useDivisor && !useByDivisor)
                 {
                         MersenneTesters = new ThreadLocal<MersenneNumberTester>(() =>
 			{
@@ -461,7 +478,7 @@ internal static class Program
 				return tester;
 			}, trackAllValues: true);
 		}
-		else
+                else if (useDivisor)
                 {
                         if (divisor == UInt128.Zero)
                         {
@@ -470,6 +487,10 @@ internal static class Program
 
                         _divisorTester = new MersenneNumberDivisorGpuTester();
                 }
+                else if (useByDivisor)
+                {
+                        _byDivisorTester = new MersenneNumberDivisorByDivisorGpuTester();
+                }
 
 		// Load RLE blacklist (optional)
 		if (!string.IsNullOrEmpty(_rleBlacklistPath))
@@ -502,6 +523,7 @@ internal static class Program
                                 kernelType,
                                 useLucas,
                                 useDivisor,
+                                useByDivisor,
                                 mersenneOnGpu,
 				useOrder,
 				useModuloWorkaround,
@@ -591,7 +613,7 @@ internal static class Program
 		// Configure batch size for GPU primality sieve
 		PrimeTester.GpuBatchSize = gpuPrimeBatch;
 
-		if (!useDivisor)
+                if (!useDivisor && !useByDivisor)
                 {
                         Console.WriteLine("Warming up orders...");
 
@@ -734,7 +756,7 @@ internal static class Program
 		Console.WriteLine("  --increment=bit|add    exponent increment method");
 		Console.WriteLine("  --threads=<value>      number of worker threads");
 		Console.WriteLine("  --block-size=<value>   values processed per thread batch");
-		Console.WriteLine("  --mersenne=pow2mod|incremental|lucas|residue|divisor  Mersenne test method");
+                Console.WriteLine("  --mersenne=pow2mod|incremental|lucas|residue|divisor|bydivisor  Mersenne test method");
 		Console.WriteLine("  --divisor=<value>     optional divisor for --mersenne=divisor mode");
 		Console.WriteLine("  --residue-max-k=<value>  max k for residue Mersenne test (q = 2*p*k + 1)");
 		Console.WriteLine("  --mersenne-device=cpu|gpu  Device for Mersenne method (default gpu)");
@@ -769,10 +791,16 @@ internal static class Program
 		Console.WriteLine("  --help, -help, --?, -?, /?   show this help message");
 	}
 
-	private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice)
+        private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice)
 	{
 		string inc = bitInc ? "bit" : "add";
-		string mers = useDivisorFlag ? "divisor" : (useLucasFlag ? "lucas" : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental"));
+                string mers = useDivisorFlag
+                                ? "divisor"
+                                : (useLucasFlag
+                                                ? "lucas"
+                                                : (useByDivisorFlag
+                                                                ? "bydivisor"
+                                                                : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental")));
 		string ntt = nttBackend == NttBackend.Staged ? "staged" : "reference";
 		string red = reduction switch { ModReductionMode.Mont64 => "mont64", ModReductionMode.Barrett128 => "barrett128", ModReductionMode.GpuUInt128 => "uint128", _ => "auto" };
 		string order = useOrder ? "order-on" : "order-off";
@@ -1113,6 +1141,11 @@ internal static class Program
 		}
 
                 searchedMersenne = true;
+                if (_useByDivisorMode)
+                {
+                        return _byDivisorTester!.IsPrime(p, divisorCyclesSearchLimit, out detailedCheck);
+                }
+
                 if (_useDivisor)
                 {
                         return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck);

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -60,5 +60,20 @@ public class MersenneNumberDivisorGpuTesterTests
         tester.IsPrime(3UL, 7UL, ulong.MaxValue, out bool exhausted).Should().BeFalse();
         exhausted.Should().BeTrue();
     }
+
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_tester_tracks_divisors_across_primes()
+    {
+        var tester = new MersenneNumberDivisorByDivisorGpuTester();
+        tester.IsPrime(5UL, 64UL, out bool exhausted).Should().BeTrue();
+        exhausted.Should().BeTrue();
+
+        tester.IsPrime(7UL, 96UL, out exhausted).Should().BeTrue();
+        exhausted.Should().BeTrue();
+
+        tester.IsPrime(11UL, 128UL, out exhausted).Should().BeFalse();
+        exhausted.Should().BeTrue();
+    }
 }
 

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -0,0 +1,307 @@
+using System;
+using System.Buffers;
+using System.Collections.Concurrent;
+using System.Collections.Generic;
+using System.Runtime.InteropServices;
+using ILGPU;
+using ILGPU.Algorithms;
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public sealed class MersenneNumberDivisorByDivisorGpuTester
+{
+        private readonly List<ulong> _divisors = new();
+        private readonly List<ulong> _remainders = new();
+        private readonly List<ulong> _lastPrimes = new();
+        private readonly object _sync = new();
+        private ulong _divisorLimit;
+
+        private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>> _updateKernelCache = new();
+        private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>> _initialKernelCache = new();
+
+        private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> GetUpdateKernel(Accelerator accelerator) =>
+                        _updateKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(UpdateKernel));
+
+        private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> GetInitialKernel(Accelerator accelerator) =>
+                        _initialKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(InitialKernel));
+
+        public bool IsPrime(ulong prime, ulong divisorLimit, out bool divisorsExhausted)
+        {
+                lock (_sync)
+                {
+                        _divisorLimit = divisorLimit;
+                        ulong allowedMax = ComputeAllowedMaxDivisor(prime, divisorLimit);
+                        if (allowedMax < 3UL)
+                        {
+                                divisorsExhausted = true;
+                                return true;
+                        }
+
+                        bool composite = UpdateExistingDivisors(prime, allowedMax);
+                        if (!composite)
+                        {
+                                composite = ExtendDivisors(prime, allowedMax);
+                        }
+
+                        divisorsExhausted = true;
+                        return !composite;
+                }
+        }
+
+        private bool UpdateExistingDivisors(ulong prime, ulong allowedMax)
+        {
+                int count = _divisors.Count;
+                if (count == 0)
+                {
+                        return false;
+                }
+
+                var gpuLease = GpuContextPool.RentPreferred(preferCpu: false);
+                var accelerator = gpuLease.Accelerator;
+                var kernel = GetUpdateKernel(accelerator);
+                var divisorsBuffer = accelerator.Allocate1D<ulong>(count);
+                var remaindersBuffer = accelerator.Allocate1D<ulong>(count);
+                var lastPrimesBuffer = accelerator.Allocate1D<ulong>(count);
+                var hitsBuffer = accelerator.Allocate1D<byte>(count);
+
+                Span<ulong> divisorsSpan = CollectionsMarshal.AsSpan(_divisors);
+                Span<ulong> remaindersSpan = CollectionsMarshal.AsSpan(_remainders);
+                Span<ulong> lastPrimesSpan = CollectionsMarshal.AsSpan(_lastPrimes);
+
+                divisorsBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(divisorsSpan), count);
+                remaindersBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(remaindersSpan), count);
+                lastPrimesBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(lastPrimesSpan), count);
+                hitsBuffer.MemSetToZero();
+
+                kernel(count, prime, divisorsBuffer.View, remaindersBuffer.View, lastPrimesBuffer.View, hitsBuffer.View);
+                accelerator.Synchronize();
+
+                remaindersBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(remaindersSpan), count);
+                lastPrimesBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(lastPrimesSpan), count);
+
+                byte[] hits = ArrayPool<byte>.Shared.Rent(count);
+                Span<byte> hitsSpan = hits.AsSpan(0, count);
+                hitsBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), count);
+
+                bool composite = false;
+                for (int i = 0; i < count; i++)
+                {
+                        if (divisorsSpan[i] > allowedMax)
+                        {
+                                break;
+                        }
+
+                        if (hitsSpan[i] != 0)
+                        {
+                                composite = true;
+                                break;
+                        }
+                }
+
+                ArrayPool<byte>.Shared.Return(hits, clearArray: true);
+
+                hitsBuffer.Dispose();
+                lastPrimesBuffer.Dispose();
+                remaindersBuffer.Dispose();
+                divisorsBuffer.Dispose();
+                gpuLease.Dispose();
+
+                return composite;
+        }
+
+        private bool ExtendDivisors(ulong prime, ulong allowedMax)
+        {
+                ulong currentMax = _divisors.Count == 0 ? 1UL : _divisors[^1];
+                ulong start = currentMax < 3UL ? 3UL : currentMax + 2UL;
+                if (start > allowedMax)
+                {
+                        return false;
+                }
+
+                List<ulong> newDivisors = new();
+                for (ulong divisor = start; divisor <= allowedMax; divisor += 2UL)
+                {
+                        newDivisors.Add(divisor);
+                }
+
+                int count = newDivisors.Count;
+                if (count == 0)
+                {
+                        return false;
+                }
+
+                var gpuLease = GpuContextPool.RentPreferred(preferCpu: false);
+                var accelerator = gpuLease.Accelerator;
+                var kernel = GetInitialKernel(accelerator);
+                var divisorsBuffer = accelerator.Allocate1D<ulong>(count);
+                var remaindersBuffer = accelerator.Allocate1D<ulong>(count);
+                var hitsBuffer = accelerator.Allocate1D<byte>(count);
+
+                Span<ulong> newDivisorsSpan = CollectionsMarshal.AsSpan(newDivisors);
+                divisorsBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(newDivisorsSpan), count);
+                hitsBuffer.MemSetToZero();
+
+                kernel(count, prime, divisorsBuffer.View, remaindersBuffer.View, hitsBuffer.View);
+                accelerator.Synchronize();
+
+                ulong[] newRemainders = ArrayPool<ulong>.Shared.Rent(count);
+                byte[] hits = ArrayPool<byte>.Shared.Rent(count);
+                Span<ulong> newRemaindersSpan = newRemainders.AsSpan(0, count);
+                Span<byte> hitsSpan = hits.AsSpan(0, count);
+                remaindersBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(newRemaindersSpan), count);
+                hitsBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), count);
+
+                bool composite = false;
+                for (int i = 0; i < count; i++)
+                {
+                        ulong divisor = newDivisorsSpan[i];
+                        ulong remainder = newRemaindersSpan[i];
+                        _divisors.Add(divisor);
+                        _remainders.Add(remainder);
+                        _lastPrimes.Add(prime);
+                        if (!composite && hitsSpan[i] != 0)
+                        {
+                                composite = true;
+                        }
+                }
+
+                ArrayPool<ulong>.Shared.Return(newRemainders, clearArray: true);
+                ArrayPool<byte>.Shared.Return(hits, clearArray: true);
+
+                hitsBuffer.Dispose();
+                remaindersBuffer.Dispose();
+                divisorsBuffer.Dispose();
+                gpuLease.Dispose();
+
+                return composite;
+        }
+
+        private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
+        {
+                if (prime <= 1UL)
+                {
+                        return 0UL;
+                }
+
+                if (prime - 1UL >= 64UL)
+                {
+                        return divisorLimit;
+                }
+
+                ulong maxByPrime = (1UL << (int)(prime - 1UL)) - 1UL;
+                return Math.Min(maxByPrime, divisorLimit);
+        }
+
+        private static void UpdateKernel(Index1D index, ulong prime, ArrayView<ulong> divisors, ArrayView<ulong> remainders, ArrayView<ulong> lastPrimes, ArrayView<byte> hits)
+        {
+                ulong divisor = divisors[index];
+                if (divisor <= 1UL || (divisor & 1UL) == 0UL)
+                {
+                        hits[index] = 0;
+                        return;
+                }
+
+                ulong lastPrime = lastPrimes[index];
+                ulong remainder = remainders[index] % divisor;
+                ulong step = prime >= lastPrime ? prime - lastPrime : 0UL;
+                if (step == 0UL)
+                {
+                        hits[index] = remainder == 0UL ? (byte)1 : (byte)0;
+                        return;
+                }
+
+                ulong pow = Pow2Mod(step, divisor);
+                ulong temp = remainder + 1UL;
+                if (temp >= divisor)
+                {
+                        temp -= divisor;
+                }
+
+                ulong product = MulMod(temp, pow, divisor);
+                ulong newRemainder = product == 0UL ? divisor - 1UL : product - 1UL;
+                bool isZero = newRemainder == 0UL;
+                remainders[index] = newRemainder;
+                lastPrimes[index] = prime;
+                hits[index] = isZero ? (byte)1 : (byte)0;
+        }
+
+        private static void InitialKernel(Index1D index, ulong prime, ArrayView<ulong> divisors, ArrayView<ulong> remainders, ArrayView<byte> hits)
+        {
+                ulong divisor = divisors[index];
+                if (divisor <= 1UL || (divisor & 1UL) == 0UL)
+                {
+                        remainders[index] = 0UL;
+                        hits[index] = 0;
+                        return;
+                }
+
+                ulong pow = Pow2Mod(prime, divisor);
+                ulong remainder = pow == 0UL ? divisor - 1UL : (pow == 1UL ? 0UL : pow - 1UL);
+                remainders[index] = remainder;
+                hits[index] = remainder == 0UL ? (byte)1 : (byte)0;
+        }
+
+        private static ulong Pow2Mod(ulong exponent, ulong modulus)
+        {
+                if (modulus <= 1UL)
+                {
+                        return 0UL;
+                }
+
+                ulong result = 1UL % modulus;
+                ulong baseVal = 2UL % modulus;
+                ulong exp = exponent;
+                while (exp > 0UL)
+                {
+                        if ((exp & 1UL) != 0UL)
+                        {
+                                result = MulMod(result, baseVal, modulus);
+                        }
+
+                        exp >>= 1;
+                        if (exp == 0UL)
+                        {
+                                break;
+                        }
+
+                        baseVal = MulMod(baseVal, baseVal, modulus);
+                }
+
+                return result;
+        }
+
+        private static ulong MulMod(ulong a, ulong b, ulong modulus)
+        {
+                if (modulus == 0UL)
+                {
+                        return 0UL;
+                }
+
+                ulong result = 0UL;
+                ulong x = a % modulus;
+                ulong y = b;
+                while (y > 0UL)
+                {
+                        if ((y & 1UL) != 0UL)
+                        {
+                                result += x;
+                                if (result >= modulus)
+                                {
+                                        result -= modulus;
+                                }
+                        }
+
+                        x <<= 1;
+                        if (x >= modulus)
+                        {
+                                x -= modulus;
+                        }
+
+                        y >>= 1;
+                }
+
+                return result;
+        }
+}
+