Merge pull request #9 from mlemanczyk/codex/implement-mersenne-divisor-check-mode

Add GPU divisor mode for Mersenne checks

EvenPerfectBitScanner.Tests/ResultsFileNameTests.cs

@@ -19,7 +19,7 @@ public class ResultsFileNameTests
             types: new[]
             {
                 typeof(bool), typeof(int), typeof(int), typeof(global::PerfectNumbers.Core.GpuKernelType),
-                typeof(bool), typeof(bool), typeof(bool), typeof(bool), typeof(bool), typeof(NttBackend), typeof(int), typeof(int), typeof(int), typeof(ulong), typeof(ModReductionMode), typeof(string), typeof(string), typeof(string)
+                typeof(bool), typeof(bool), typeof(bool), typeof(bool), typeof(bool), typeof(bool), typeof(NttBackend), typeof(int), typeof(int), typeof(int), typeof(ulong), typeof(ModReductionMode), typeof(string), typeof(string), typeof(string)
             },
             modifiers: null)!;
 
@@ -30,6 +30,7 @@ public class ResultsFileNameTests
             1,               // block
             global::PerfectNumbers.Core.GpuKernelType.Incremental,
             false,           // useLucasFlag
+            false,           // useDivisorFlag
             true,            // mersenneOnGpu (unused inside name content, devices passed explicitly below)
             false,           // useOrder
             false,           // useModWorkaround

EvenPerfectBitScanner/Program.cs

@@ -26,6 +26,9 @@ internal static class Program
 	private static bool _primeFoundAfterInit;
 	private static bool _useGcdFilter;
         private static bool _useResidue;
+        private static bool _useDivisor;
+        private static ulong _divisor;
+        private static MersenneNumberDivisorGpuTester? _divisorTester;
 	private static ulong? _orderWarmupLimitOverride;
 	private static unsafe delegate*<ulong, ref ulong, ulong> _transformP;
 	private static long _state;
@@ -58,6 +61,8 @@ internal static class Program
 		bool useBitTransform = false;
                 bool useLucas = false;
                 bool useResidue = false;    // M_p test via residue divisors (replaces LL/incremental)
+                bool useDivisor = false;     // M_p divisibility by specific divisor
+                ulong divisor = 0UL;
 									// Device routing
 		bool useGpuCycles = true;
 		bool mersenneOnGpu = true;   // controls Lucas/incremental/pow2mod device
@@ -113,11 +118,22 @@ internal static class Program
                                         useResidue = true;
                                         useLucas = false;
                                 }
+                                else if (value.Equals("divisor", StringComparison.OrdinalIgnoreCase))
+                                {
+                                        useDivisor = true;
+                                        useLucas = false;
+                                        useResidue = false;
+                                }
                                 else
                                 {
                                         kernelType = GpuKernelType.Incremental;
                                 }
                         }
+                        else if (arg.StartsWith("--divisor=", StringComparison.OrdinalIgnoreCase))
+                        {
+                                int eq = arg.IndexOf('=');
+                                divisor = ulong.Parse(arg.AsSpan(eq + 1));
+                        }
 			else if (arg.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase))
 			{
 				int eq = arg.IndexOf('=');
@@ -341,6 +357,12 @@ internal static class Program
                         return;
                 }
 
+                if (useDivisor && divisor == 0UL)
+                {
+                        Console.WriteLine("Divisor mode requires --divisor=<value>.");
+                        return;
+                }
+
 		// Apply GPU prime sieve runtime configuration
 		GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads);
 		PerfectNumbers.Core.PrimeTester.GpuBatchSize = Math.Max(1, gpuPrimeBatch);
@@ -389,11 +411,15 @@ internal static class Program
                 Console.WriteLine("Divisor cycles are ready");
 
                 _useResidue = useResidue;
+                _useDivisor = useDivisor;
+                _divisor = divisor;
                 _transformP = useBitTransform ? &TransformPBit : &TransformPAdd;
                 PResidue = new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, initialNumber: currentP, initialized: true), trackAllValues: true);
                 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)
+                {
                         MersenneTesters = new ThreadLocal<MersenneNumberTester>(() =>
                         {
                                 var tester = new MersenneNumberTester(
@@ -414,6 +440,11 @@ internal static class Program
 
                                 return tester;
                         }, trackAllValues: true);
+                }
+                else
+                {
+                        _divisorTester = new MersenneNumberDivisorGpuTester();
+                }
 
 		// Load RLE blacklist (optional)
 		if (!string.IsNullOrEmpty(_rleBlacklistPath))
@@ -445,6 +476,7 @@ internal static class Program
                         blockSize,
                         kernelType,
                         useLucas,
+                        useDivisor,
                         mersenneOnGpu,
                         useOrder,
                         useModuloWorkaround,
@@ -523,6 +555,7 @@ internal static class Program
                         blockSize,
                         kernelType,
                         useLucas,
+                        useDivisor,
                         mersenneOnGpu,
                         useOrder,
                         useModuloWorkaround,
@@ -537,6 +570,8 @@ internal static class Program
                         GpuContextPool.ForceCpu ? "cpu" : "gpu",
                         orderOnGpu ? "gpu" : "cpu");
 
+                if (!useDivisor)
+                {
                         Console.WriteLine("Warming up orders...");
 
                         threadCount = Math.Max(1, threadCount);
@@ -545,6 +580,7 @@ internal static class Program
                         {
                                 MersenneTesters.Value!.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL);
                         }
+                }
 
 		Console.WriteLine("Starting scan...");
 		_state = ((long)currentP << 3) | (long)remainder;
@@ -653,7 +689,8 @@ 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  Mersenne test method");
+                Console.WriteLine("  --mersenne=pow2mod|incremental|lucas|residue|divisor  Mersenne test method");
+                Console.WriteLine("  --divisor=<value>     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)");
                 Console.WriteLine("  --primes-device=cpu|gpu    Device for prime-scan kernels (default gpu)");
@@ -672,6 +709,10 @@ internal static class Program
 		Console.WriteLine("  --zero-conj=<f>:<r>     hard reject if zero_fraction(p) > f AND max_zero_block >= r (default off)");
                 Console.WriteLine("  --results-dir=<path>   directory for results file");
                 Console.WriteLine("  --results-prefix=<text> prefix to prepend to results filename");
+                Console.WriteLine("  --divisor-cycles=<path>       divisor cycles data file");
+                Console.WriteLine("  --divisor-cycles-device=cpu|gpu  device for cycles generation (default gpu)");
+                Console.WriteLine("  --divisor-cycles-batch-size=<value> batch size for cycles generation (default 512)");
+                Console.WriteLine("  --divisor-cycles-continue  continue divisor cycles generation");
 
 		Console.WriteLine("  --use-order            test primality via q order");
 		Console.WriteLine("  --workaround-mod       avoid '%' operator on the GPU");
@@ -682,10 +723,10 @@ 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 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 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 = useLucasFlag ? "lucas" : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental");
+                string mers = useDivisorFlag ? "divisor" : (useLucasFlag ? "lucas" : (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";
@@ -929,6 +970,12 @@ internal static class Program
 		}
 
                 searchedMersenne = true;
+                if (_useDivisor)
+                {
+                        detailedCheck = _divisorTester!.IsDivisible(p, _divisor);
+                        return detailedCheck;
+                }
+
                 detailedCheck = MersenneTesters.Value!.IsMersennePrime(p);
                 return detailedCheck;
         }

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -0,0 +1,20 @@
+using FluentAssertions;
+using PerfectNumbers.Core.Gpu;
+using Xunit;
+
+namespace PerfectNumbers.Core.Tests;
+
+public class MersenneNumberDivisorGpuTesterTests
+{
+    [Theory]
+    [InlineData(3UL, 7UL, true)]
+    [InlineData(11UL, 23UL, true)]
+    [InlineData(7UL, 35UL, false)]
+    [InlineData(13UL, 23UL, false)]
+    [Trait("Category", "Fast")]
+    public void IsDivisible_returns_expected(ulong exponent, ulong divisor, bool expected)
+    {
+        var tester = new MersenneNumberDivisorGpuTester();
+        tester.IsDivisible(exponent, divisor).Should().Be(expected);
+    }
+}

PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs

@@ -0,0 +1,57 @@
+using System.Collections.Concurrent;
+using ILGPU;
+using ILGPU.Algorithms;
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public sealed class MersenneNumberDivisorGpuTester
+{
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ulong, ArrayView<byte>>> _kernelCache = new();
+
+    private Action<Index1D, ulong, ulong, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
+        _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<byte>>(Kernel));
+
+    public bool IsDivisible(ulong exponent, ulong divisor)
+    {
+        var gpu = GpuContextPool.RentPreferred(preferCpu: false);
+        var accelerator = gpu.Accelerator;
+        var kernel = GetKernel(accelerator);
+        using var resultBuffer = accelerator.Allocate1D<byte>(1);
+        kernel(1, exponent, divisor, resultBuffer.View);
+        accelerator.Synchronize();
+        bool divisible = resultBuffer.GetAsArray1D()[0] != 0;
+        gpu.Dispose();
+        return divisible;
+    }
+
+    private static void Kernel(Index1D _, ulong exponent, ulong divisor, ArrayView<byte> result)
+    {
+        if (divisor <= 1UL)
+        {
+            result[0] = 0;
+            return;
+        }
+
+        int x = 63 - XMath.LeadingZeroCount(divisor);
+        if (x == 0)
+        {
+            var mod = new GpuUInt128(0UL, divisor);
+            var temp = GpuUInt128.Pow2Mod(exponent, mod);
+            result[0] = temp.High == 0UL && temp.Low == 1UL ? (byte)1 : (byte)0;
+            return;
+        }
+
+        ulong ux = (ulong)x;
+        ulong q = exponent / ux;
+        ulong r = exponent % ux;
+        ulong pow2x = 1UL << x;
+        ulong y = divisor - pow2x;
+        var modVal = new GpuUInt128(0UL, divisor);
+        var baseVal = new GpuUInt128(0UL, pow2x); // 2^x ≡ -y (mod divisor)
+        var part1 = baseVal.ModPow(q, modVal);
+        var part2 = GpuUInt128.Pow2Mod(r, modVal);
+        var rem = part1.MulMod(part2, modVal);
+        result[0] = rem.High == 0UL && rem.Low == 1UL ? (byte)1 : (byte)0;
+    }
+}