diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index d6c3d50..85f5bab 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -28,10 +28,12 @@ internal static class Program private static int _primeCount; private static bool _primeFoundAfterInit; private static bool _useGcdFilter; - private static bool _useDivisor; - private static bool _useResidueMode; - private static UInt128 _divisor; - private static MersenneNumberDivisorGpuTester? _divisorTester; + 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* _transformP; private static long _state; @@ -66,8 +68,9 @@ internal static class Program bool showHelp = false; 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 + 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; @@ -125,12 +128,19 @@ internal static class Program useResidue = true; useLucas = false; } - else if (value.Equals("divisor", StringComparison.OrdinalIgnoreCase)) - { - useDivisor = true; - useLucas = false; - useResidue = false; - } + else if (value.Equals("divisor", StringComparison.OrdinalIgnoreCase)) + { + useDivisor = true; + useLucas = false; + useResidue = false; + } + else if (value.Equals("bydivisor", StringComparison.OrdinalIgnoreCase)) + { + useByDivisor = true; + useLucas = false; + useResidue = false; + useDivisor = false; + } else { kernelType = GpuKernelType.Incremental; @@ -364,13 +374,19 @@ internal static class Program { continueCyclesGeneration = true; } - } + } - if (showHelp) - { - PrintHelp(); - return; - } + if (useByDivisor && string.IsNullOrEmpty(filterFile)) + { + Console.WriteLine("--mersenne=bydivisor requires --filter-p=."); + return; + } + + if (showHelp) + { + PrintHelp(); + return; + } // Apply GPU prime sieve runtime configuration GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); @@ -419,9 +435,16 @@ internal static class Program Console.WriteLine("Divisor cycles are ready"); - _useDivisor = useDivisor; - _useResidueMode = useResidue; - _divisor = divisor; + if (useByDivisor) + { + threadCount = 1; + blockSize = 1; + } + + _useDivisor = useDivisor; + _useResidueMode = useResidue; + _useByDivisorMode = useByDivisor; + _divisor = divisor; if (useBitTransform) { _transformP = &TransformPBit; @@ -438,9 +461,9 @@ internal static class Program PrimeTesters = new ThreadLocal(() => 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(() => + if (!useDivisor && !useByDivisor) + { + MersenneTesters = new ThreadLocal(() => { var tester = new MersenneNumberTester( useIncremental: !useLucas, @@ -461,15 +484,19 @@ internal static class Program return tester; }, trackAllValues: true); } - else - { - if (divisor == UInt128.Zero) - { - MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); - } + else if (useDivisor) + { + if (divisor == UInt128.Zero) + { + MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); + } - _divisorTester = new MersenneNumberDivisorGpuTester(); - } + _divisorTester = new MersenneNumberDivisorGpuTester(); + } + else if (useByDivisor) + { + _byDivisorTester = new MersenneNumberDivisorByDivisorGpuTester(); + } // Load RLE blacklist (optional) if (!string.IsNullOrEmpty(_rleBlacklistPath)) @@ -495,14 +522,15 @@ internal static class Program Console.WriteLine("Initialization..."); // Compose a results file name that encodes configuration (before opening file) - var builtName = BuildResultsFileName( - useBitTransform, - threadCount, - blockSize, - kernelType, - useLucas, - useDivisor, - mersenneOnGpu, + var builtName = BuildResultsFileName( + useBitTransform, + threadCount, + blockSize, + kernelType, + useLucas, + useDivisor, + useByDivisor, + mersenneOnGpu, useOrder, useModuloWorkaround, _useGcdFilter, @@ -544,12 +572,12 @@ internal static class Program File.WriteAllText(ResultsFileName, $"p,searchedMersenne,detailedCheck,passedAllTests{Environment.NewLine}"); } - bool useFilter = !string.IsNullOrEmpty(filterFile); - HashSet filter = []; - ulong maxP = 0UL; - if (useFilter) - { - Console.WriteLine("Loading filter..."); + bool useFilter = !string.IsNullOrEmpty(filterFile); + HashSet filter = []; + ulong maxP = 0UL; + if (useFilter) + { + Console.WriteLine("Loading filter..."); ulong[] localFilter = new ulong[1024]; int count = 0; LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) => @@ -574,15 +602,26 @@ internal static class Program if (count > 0) { filter.AddRange(localFilter[..count]); - } - // Restore this if you want to use List instead of HashSet - // filter.Sort(); - } + } + // Restore this if you want to use List instead of HashSet + // filter.Sort(); + } - if (_primeCount >= 2) - { - _primeFoundAfterInit = true; - } + 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; + } _outputBuilder = StringBuilderPool.Rent(); @@ -591,11 +630,11 @@ internal static class Program // Configure batch size for GPU primality sieve PrimeTester.GpuBatchSize = gpuPrimeBatch; - if (!useDivisor) - { - Console.WriteLine("Warming up orders..."); + if (!useDivisor && !useByDivisor) + { + Console.WriteLine("Warming up orders..."); - threadCount = Math.Max(1, threadCount); + threadCount = Math.Max(1, threadCount); _ = MersenneTesters.Value; if (!useLucas) { @@ -734,7 +773,7 @@ internal static class Program Console.WriteLine(" --increment=bit|add exponent increment method"); Console.WriteLine(" --threads= number of worker threads"); Console.WriteLine(" --block-size= 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= optional divisor for --mersenne=divisor mode"); Console.WriteLine(" --residue-max-k= 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= process only p from previous run results"); + Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor)"); Console.WriteLine(" --write-batch-size= 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"; @@ -1112,11 +1157,16 @@ internal static class Program return false; } - searchedMersenne = true; - if (_useDivisor) - { - return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); - } + searchedMersenne = true; + if (_useByDivisorMode) + { + return _byDivisorTester!.IsPrime(p, out detailedCheck); + } + + if (_useDivisor) + { + return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); + } detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); return detailedCheck; diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 18f1258..d50edf3 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/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(); + } + + [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(); + } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs new file mode 100644 index 0000000..6133440 --- /dev/null +++ b/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 _divisors = new(); + private readonly List _remainders = new(); + private readonly List _lastPrimes = new(); + private readonly object _sync = new(); + private ulong _divisorLimit; + private bool _isConfigured; + private ulong _lastStatusDivisor; + + private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView>> _updateKernelCache = new(); + private readonly ConcurrentDictionary, ArrayView, ArrayView>> _initialKernelCache = new(); + + private Action, ArrayView, ArrayView, ArrayView> GetUpdateKernel(Accelerator accelerator) => + _updateKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(UpdateKernel)); + + private Action, ArrayView, ArrayView> GetInitialKernel(Accelerator accelerator) => + _initialKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(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(count); + var remaindersBuffer = accelerator.Allocate1D(count); + var lastPrimesBuffer = accelerator.Allocate1D(count); + var hitsBuffer = accelerator.Allocate1D(count); + + Span divisorsSpan = CollectionsMarshal.AsSpan(_divisors); + Span remaindersSpan = CollectionsMarshal.AsSpan(_remainders); + Span 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.Shared.Rent(count); + Span 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.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 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(count); + var remaindersBuffer = accelerator.Allocate1D(count); + var hitsBuffer = accelerator.Allocate1D(count); + + Span 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.Shared.Rent(count); + byte[] hits = ArrayPool.Shared.Rent(count); + Span newRemaindersSpan = newRemainders.AsSpan(0, count); + Span 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.Shared.Return(newRemainders, clearArray: true); + ArrayPool.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 divisors, ArrayView remainders, ArrayView lastPrimes, ArrayView 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 divisors, ArrayView remainders, ArrayView 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; + } +} +