From 596bb7e38fff63e61cde91c434a18a2d9aaec62b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 10 Sep 2025 22:03:17 +0200 Subject: [PATCH] Add GPU divisor tester --- .../ResultsFileNameTests.cs | 3 +- EvenPerfectBitScanner/Program.cs | 271 ++++++++++-------- .../MersenneNumberDivisorGpuTesterTests.cs | 20 ++ .../Gpu/MersenneNumberDivisorGpuTester.cs | 57 ++++ 4 files changed, 238 insertions(+), 113 deletions(-) create mode 100644 PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs create mode 100644 PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs diff --git a/EvenPerfectBitScanner.Tests/ResultsFileNameTests.cs b/EvenPerfectBitScanner.Tests/ResultsFileNameTests.cs index 806ebbb..1b4c5c7 100644 --- a/EvenPerfectBitScanner.Tests/ResultsFileNameTests.cs +++ b/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 diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index dae6d58..05f0d2a 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -25,7 +25,10 @@ internal static class Program private static int _primeCount; private static bool _primeFoundAfterInit; private static bool _useGcdFilter; - private static bool _useResidue; + private static bool _useResidue; + private static bool _useDivisor; + private static ulong _divisor; + private static MersenneNumberDivisorGpuTester? _divisorTester; private static ulong? _orderWarmupLimitOverride; private static unsafe delegate* _transformP; private static long _state; @@ -56,8 +59,10 @@ internal static class Program bool useOrder = false; bool showHelp = false; bool useBitTransform = false; - bool useLucas = false; - bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) + 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 @@ -100,24 +105,35 @@ internal static class Program { int eq = arg.IndexOf('='); ReadOnlySpan value = arg.AsSpan(eq + 1); - if (value.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) - { - kernelType = GpuKernelType.Pow2Mod; - } - else if (value.Equals("lucas", StringComparison.OrdinalIgnoreCase)) - { - useLucas = true; - } - else if (value.Equals("residue", StringComparison.OrdinalIgnoreCase)) - { - useResidue = true; - useLucas = false; - } - else - { - kernelType = GpuKernelType.Incremental; - } - } + if (value.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) + { + kernelType = GpuKernelType.Pow2Mod; + } + else if (value.Equals("lucas", StringComparison.OrdinalIgnoreCase)) + { + useLucas = true; + } + else if (value.Equals("residue", StringComparison.OrdinalIgnoreCase)) + { + 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('='); @@ -335,11 +351,17 @@ internal static class Program } } - if (showHelp) - { - PrintHelp(); - return; - } + if (showHelp) + { + PrintHelp(); + return; + } + + if (useDivisor && divisor == 0UL) + { + Console.WriteLine("Divisor mode requires --divisor=."); + return; + } // Apply GPU prime sieve runtime configuration GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); @@ -386,34 +408,43 @@ internal static class Program MersenneDivisorCycles.Shared.LoadFrom(cyclesPath); } - Console.WriteLine("Divisor cycles are ready"); + Console.WriteLine("Divisor cycles are ready"); - _useResidue = useResidue; - _transformP = useBitTransform ? &TransformPBit : &TransformPAdd; - PResidue = new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, initialNumber: currentP, initialized: true), trackAllValues: true); - 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 - MersenneTesters = new ThreadLocal(() => - { - var tester = new MersenneNumberTester( - useIncremental: !useLucas, - useOrderCache: false, - kernelType: kernelType, - useModuloWorkaround: useModuloWorkaround, - useOrder: useOrder, - useGpuLucas: mersenneOnGpu, - useGpuScan: mersenneOnGpu, - useGpuOrder: orderOnGpu, - useResidue: useResidue, - maxK: residueKMax); - if (!useLucas) - { - tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); - } + _useResidue = useResidue; + _useDivisor = useDivisor; + _divisor = divisor; + _transformP = useBitTransform ? &TransformPBit : &TransformPAdd; + PResidue = new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, initialNumber: currentP, initialized: true), trackAllValues: true); + 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(() => + { + var tester = new MersenneNumberTester( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useModuloWorkaround: useModuloWorkaround, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu, + useResidue: useResidue, + maxK: residueKMax); + if (!useLucas) + { + tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); + } - return tester; - }, trackAllValues: true); + return tester; + }, trackAllValues: true); + } + else + { + _divisorTester = new MersenneNumberDivisorGpuTester(); + } // Load RLE blacklist (optional) if (!string.IsNullOrEmpty(_rleBlacklistPath)) @@ -439,25 +470,26 @@ 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, - mersenneOnGpu, - useOrder, - useModuloWorkaround, - _useGcdFilter, - NttGpuMath.GpuTransformBackend, - gpuPrimeThreads, - sliceSize, - scanBatchSize, - _orderWarmupLimitOverride ?? 5_000_000UL, - NttGpuMath.ReductionMode, - mersenneOnGpu ? "gpu" : "cpu", - (GpuContextPool.ForceCpu ? "cpu" : "gpu"), - orderOnGpu ? "gpu" : "cpu"); + var builtName = BuildResultsFileName( + useBitTransform, + threadCount, + blockSize, + kernelType, + useLucas, + useDivisor, + mersenneOnGpu, + useOrder, + useModuloWorkaround, + _useGcdFilter, + NttGpuMath.GpuTransformBackend, + gpuPrimeThreads, + sliceSize, + scanBatchSize, + _orderWarmupLimitOverride ?? 5_000_000UL, + NttGpuMath.ReductionMode, + mersenneOnGpu ? "gpu" : "cpu", + (GpuContextPool.ForceCpu ? "cpu" : "gpu"), + orderOnGpu ? "gpu" : "cpu"); if (!string.IsNullOrEmpty(_resultsPrefix)) { @@ -517,34 +549,38 @@ internal static class Program PerfectNumbers.Core.PrimeTester.GpuBatchSize = gpuPrimeBatch; // Compose a results file name that encodes configuration - ResultsFileName = BuildResultsFileName( - useBitTransform, - threadCount, - blockSize, - kernelType, - useLucas, - mersenneOnGpu, - useOrder, - useModuloWorkaround, - _useGcdFilter, - NttGpuMath.GpuTransformBackend, - gpuPrimeThreads, - sliceSize, - scanBatchSize, - _orderWarmupLimitOverride ?? 5_000_000UL, - NttGpuMath.ReductionMode, - mersenneOnGpu ? "gpu" : "cpu", - GpuContextPool.ForceCpu ? "cpu" : "gpu", - orderOnGpu ? "gpu" : "cpu"); + ResultsFileName = BuildResultsFileName( + useBitTransform, + threadCount, + blockSize, + kernelType, + useLucas, + useDivisor, + mersenneOnGpu, + useOrder, + useModuloWorkaround, + _useGcdFilter, + NttGpuMath.GpuTransformBackend, + gpuPrimeThreads, + sliceSize, + scanBatchSize, + _orderWarmupLimitOverride ?? 5_000_000UL, + NttGpuMath.ReductionMode, + mersenneOnGpu ? "gpu" : "cpu", + GpuContextPool.ForceCpu ? "cpu" : "gpu", + orderOnGpu ? "gpu" : "cpu"); - Console.WriteLine("Warming up orders..."); + if (!useDivisor) + { + Console.WriteLine("Warming up orders..."); - threadCount = Math.Max(1, threadCount); - _ = MersenneTesters.Value; - if (!useLucas) - { - MersenneTesters.Value!.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); - } + threadCount = Math.Max(1, threadCount); + _ = MersenneTesters.Value; + if (!useLucas) + { + MersenneTesters.Value!.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); + } + } Console.WriteLine("Starting scan..."); _state = ((long)currentP << 3) | (long)remainder; @@ -653,11 +689,12 @@ 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 Mersenne test method"); - 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)"); - Console.WriteLine(" --primes-device=cpu|gpu Device for prime-scan kernels (default gpu)"); - Console.WriteLine(" --gpu-prime-batch= Batch size for GPU primality sieve (default 262144)"); + Console.WriteLine(" --mersenne=pow2mod|incremental|lucas|residue|divisor Mersenne test method"); + Console.WriteLine(" --divisor= 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)"); + Console.WriteLine(" --primes-device=cpu|gpu Device for prime-scan kernels (default gpu)"); + Console.WriteLine(" --gpu-prime-batch= Batch size for GPU primality sieve (default 262144)"); Console.WriteLine(" --order-device=cpu|gpu Device for order computations (default gpu)"); Console.WriteLine(" --ntt=reference|staged GPU NTT backend (default staged)"); Console.WriteLine(" --mod-reduction=auto|uint128|mont64|barrett128 staged NTT reduction (default auto)"); @@ -670,8 +707,12 @@ internal static class Program Console.WriteLine(" --rle-only-last7=true|false apply RLE only when p % 10 == 7 (default true)"); Console.WriteLine(" --zero-hard= hard reject if zero_fraction(p) > f (default off)"); Console.WriteLine(" --zero-conj=: hard reject if zero_fraction(p) > f AND max_zero_block >= r (default off)"); - Console.WriteLine(" --results-dir= directory for results file"); - Console.WriteLine(" --results-prefix= prefix to prepend to results filename"); + Console.WriteLine(" --results-dir= directory for results file"); + Console.WriteLine(" --results-prefix= prefix to prepend to results filename"); + Console.WriteLine(" --divisor-cycles= divisor cycles data file"); + Console.WriteLine(" --divisor-cycles-device=cpu|gpu device for cycles generation (default gpu)"); + Console.WriteLine(" --divisor-cycles-batch-size= 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) - { - string inc = bitInc ? "bit" : "add"; - string mers = useLucasFlag ? "lucas" : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental"); + 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 = 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"; @@ -928,10 +969,16 @@ internal static class Program return false; } - searchedMersenne = true; - detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); - return detailedCheck; - } + searchedMersenne = true; + if (_useDivisor) + { + detailedCheck = _divisorTester!.IsDivisible(p, _divisor); + return detailedCheck; + } + + detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); + return detailedCheck; + } // Use ModResidueTracker with a small set of primes to prefilter composite p. private static bool IsCompositeByResidues(ulong p) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs new file mode 100644 index 0000000..d50a48a --- /dev/null +++ b/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); + } +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs new file mode 100644 index 0000000..9835a26 --- /dev/null +++ b/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>> _kernelCache = new(); + + private Action> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(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(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 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; + } +}