Merge pull request #54 from mlemanczyk/codex/implement-bydivisor-method-for-prime-testing

Add GPU by-divisor scanning mode

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;
@@ -366,6 +376,12 @@ internal static class Program
 			}
                 }
 
+                if (useByDivisor && string.IsNullOrEmpty(filterFile))
+                {
+                        Console.WriteLine("--mersenne=bydivisor requires --filter-p=<path>.");
+                        return;
+                }
+
                 if (showHelp)
                 {
                         PrintHelp();
@@ -419,8 +435,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 +461,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 +484,7 @@ internal static class Program
 				return tester;
 			}, trackAllValues: true);
 		}
-		else
+                else if (useDivisor)
                 {
                         if (divisor == UInt128.Zero)
                         {
@@ -470,6 +493,10 @@ internal static class Program
 
                         _divisorTester = new MersenneNumberDivisorGpuTester();
                 }
+                else if (useByDivisor)
+                {
+                        _byDivisorTester = new MersenneNumberDivisorByDivisorGpuTester();
+                }
 
 		// Load RLE blacklist (optional)
 		if (!string.IsNullOrEmpty(_rleBlacklistPath))
@@ -502,6 +529,7 @@ internal static class Program
                                 kernelType,
                                 useLucas,
                                 useDivisor,
+                                useByDivisor,
                                 mersenneOnGpu,
 				useOrder,
 				useModuloWorkaround,
@@ -579,6 +607,17 @@ internal static class Program
                         // filter.Sort();
                 }
 
+                if (useByDivisor)
+                {
+                        if (maxP <= 1UL)
+                        {
+                                Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor.");
+                                return;
+                        }
+
+                        _byDivisorTester!.ConfigureFromMaxPrime(maxP);
+                }
+
                 if (_primeCount >= 2)
                 {
                         _primeFoundAfterInit = true;
@@ -591,7 +630,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 +773,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)");
@@ -763,16 +802,22 @@ internal static class Program
 		Console.WriteLine("  --use-order            test primality via q order");
 		Console.WriteLine("  --workaround-mod       avoid '%' operator on the GPU");
 		// mod-automaton removed
-		Console.WriteLine("  --filter-p=<path>      process only p from previous run results");
+                Console.WriteLine("  --filter-p=<path>      process only p from previous run results (required for --mersenne=bydivisor)");
 		Console.WriteLine("  --write-batch-size=<value> overwrite frequency of disk writes (default 100 lines)");
 		Console.WriteLine("  --gcd-filter           enable early sieve based on GCD");
 		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 +1158,11 @@ internal static class Program
 		}
 
                 searchedMersenne = true;
+                if (_useByDivisorMode)
+                {
+                        return _byDivisorTester!.IsPrime(p, out detailedCheck);
+                }
+
                 if (_useDivisor)
                 {
                         return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck);

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -1,5 +1,6 @@
 using FluentAssertions;
 using PerfectNumbers.Core.Gpu;
+using System;
 using System.Numerics;
 using System.Reflection;
 using Xunit;
@@ -60,5 +61,46 @@ 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.ConfigureFromMaxPrime(11UL);
+
+        tester.IsPrime(5UL, out bool exhausted).Should().BeTrue();
+        exhausted.Should().BeTrue();
+
+        tester.IsPrime(7UL, out exhausted).Should().BeTrue();
+        exhausted.Should().BeTrue();
+
+        tester.IsPrime(11UL, out exhausted).Should().BeFalse();
+        exhausted.Should().BeTrue();
+    }
+
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_tester_requires_configuration()
+    {
+        var tester = new MersenneNumberDivisorByDivisorGpuTester();
+        Action act = () => tester.IsPrime(5UL, out _);
+        act.Should().Throw<InvalidOperationException>();
+    }
+
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_tester_respects_filter_based_limit()
+    {
+        var tester = new MersenneNumberDivisorByDivisorGpuTester();
+        tester.ConfigureFromMaxPrime(5UL);
+
+        tester.IsPrime(5UL, out _).Should().BeTrue();
+        tester.IsPrime(7UL, out bool exhausted).Should().BeTrue();
+        exhausted.Should().BeTrue();
+
+        tester.IsPrime(11UL, out bool divisorsExhausted).Should().BeTrue();
+        divisorsExhausted.Should().BeTrue();
+    }
 }
 

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -0,0 +1,371 @@
+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 bool _isConfigured;
+        private ulong _lastStatusDivisor;
+
+        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 void ConfigureFromMaxPrime(ulong maxPrime)
+        {
+                lock (_sync)
+                {
+                        if (_divisors.Count > 0)
+                        {
+                                _divisors.Clear();
+                                _remainders.Clear();
+                                _lastPrimes.Clear();
+                        }
+
+                        _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
+                        _lastStatusDivisor = 0UL;
+                        _isConfigured = true;
+                }
+        }
+
+        public bool IsPrime(ulong prime, out bool divisorsExhausted)
+        {
+                lock (_sync)
+                {
+                        if (!_isConfigured)
+                        {
+                                throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
+                        }
+
+                        ulong allowedMax = ComputeAllowedMaxDivisor(prime);
+                        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;
+                ulong lastProcessed = 0UL;
+                for (int i = 0; i < count; i++)
+                {
+                        ulong divisor = divisorsSpan[i];
+                        if (divisor > allowedMax)
+                        {
+                                break;
+                        }
+
+                        lastProcessed = divisor;
+                        if (hitsSpan[i] != 0)
+                        {
+                                composite = true;
+                                break;
+                        }
+                }
+
+                if (lastProcessed != 0UL)
+                {
+                        ReportStatus(lastProcessed);
+                }
+
+                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;
+                ulong lastProcessed = 0UL;
+                for (int i = 0; i < count; i++)
+                {
+                        ulong divisor = newDivisorsSpan[i];
+                        ulong remainder = newRemaindersSpan[i];
+                        _divisors.Add(divisor);
+                        _remainders.Add(remainder);
+                        _lastPrimes.Add(prime);
+                        lastProcessed = divisor;
+                        if (!composite && hitsSpan[i] != 0)
+                        {
+                                composite = true;
+                        }
+                }
+
+                if (lastProcessed != 0UL)
+                {
+                        ReportStatus(lastProcessed);
+                }
+
+                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 void ReportStatus(ulong divisor)
+        {
+                if (divisor <= _lastStatusDivisor)
+                {
+                        return;
+                }
+
+                _lastStatusDivisor = divisor;
+                Console.WriteLine($"...processed by-divisor candidate = {divisor}");
+        }
+
+        private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
+        {
+                if (maxPrime <= 1UL)
+                {
+                        return 0UL;
+                }
+
+                if (maxPrime - 1UL >= 64UL)
+                {
+                        return ulong.MaxValue;
+                }
+
+                return (1UL << (int)(maxPrime - 1UL)) - 1UL;
+        }
+
+        private ulong ComputeAllowedMaxDivisor(ulong prime)
+        {
+                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;
+        }
+}
+