From a296e62d4daac2a2693d0397249ddd16175dfab9 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 26 Nov 2025 22:03:52 +0100 Subject: [PATCH] --min-k command-line argument added. Partial factoring switched back to CPU. --- EvenPerfectBitScanner/CliArguments.cs | 25 +- EvenPerfectBitScanner/Program.cs | 1 + ...rsenneNumberDivisorByDivisorTesterTests.cs | 2 + ...MersenneNumberDivisorByDivisorCpuTester.cs | 171 +++++++-- PerfectNumbers.Core/Gpu/AcceleratorPool.cs | 20 +- .../PrimeOrderCalculatorAccelerator.cs | 25 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 351 +++++++++++++----- .../IMersenneNumberDivisorByDivisorTester.cs | 2 + PerfectNumbers.Core/PerfectNumberConstants.cs | 2 +- .../PrimeOrderCalculator.Cpu.cs | 46 +-- 10 files changed, 485 insertions(+), 160 deletions(-) diff --git a/EvenPerfectBitScanner/CliArguments.cs b/EvenPerfectBitScanner/CliArguments.cs index 8b64344..0fb64ec 100644 --- a/EvenPerfectBitScanner/CliArguments.cs +++ b/EvenPerfectBitScanner/CliArguments.cs @@ -50,6 +50,7 @@ internal readonly struct CliArguments internal readonly double ZeroFractionHard; internal readonly double ZeroFractionConjecture; internal readonly int MaxZeroConjecture; + internal readonly ulong MinK; private CliArguments( ulong startPrime, @@ -93,7 +94,8 @@ internal readonly struct CliArguments bool rleOnlyLast7, double zeroFractionHard, double zeroFractionConjecture, - int maxZeroConjecture) + int maxZeroConjecture, + ulong minK) { StartPrime = startPrime; StartPrimeProvided = startPrimeProvided; @@ -137,6 +139,7 @@ internal readonly struct CliArguments ZeroFractionHard = zeroFractionHard; ZeroFractionConjecture = zeroFractionConjecture; MaxZeroConjecture = maxZeroConjecture; + MinK = minK; } internal bool UseFilter => !TestMode && !string.IsNullOrEmpty(FilterFile); @@ -185,6 +188,7 @@ internal readonly struct CliArguments double zeroFractionHard = -1.0; double zeroFractionConjecture = -1.0; int maxZeroConjecture = -1; + ulong minK = 1UL; foreach (string argument in args) { @@ -298,6 +302,21 @@ internal readonly struct CliArguments continue; } + if (argument.StartsWith("--min-k=", StringComparison.OrdinalIgnoreCase)) + { + if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--min-k=".Length), out ulong parsedMinK) && parsedMinK > 0UL) + { + minK = parsedMinK; + } + else + { + Console.WriteLine("Invalid value for --min-k."); + return default; + } + + continue; + } + if (argument.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase)) { useMersenneOnGpu = !argument.AsSpan("--mersenne-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); @@ -568,7 +587,8 @@ internal readonly struct CliArguments rleOnlyLast7, zeroFractionHard, zeroFractionConjecture, - maxZeroConjecture); + maxZeroConjecture, + minK); } internal static void PrintHelp() @@ -609,6 +629,7 @@ internal readonly struct CliArguments 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(" --min-k= starting k for --mersenne=bydivisor divisor scan (default 1)"); Console.WriteLine(" --help, -help, --?, -?, /? show this help message"); } } diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 689a820..559c6cd 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -217,6 +217,7 @@ internal static class Program _byDivisorTester = mersenneOnGpu ? new MersenneNumberDivisorByDivisorGpuTester() : new MersenneNumberDivisorByDivisorCpuTester(); + _byDivisorTester.MinK = _cliArguments.MinK; _byDivisorTester.BatchSize = scanBatchSize; } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs index fe78197..af1e0cc 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs @@ -46,6 +46,8 @@ public class MersenneNumberDivisorByDivisorTesterTests public int BatchSize { get; set; } + public ulong MinK { get; set; } + public ulong DivisorLimit => ConfiguredMaxPrime + 100UL; public void ConfigureFromMaxPrime(ulong maxPrime) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index db6db69..2d0832f 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -8,6 +8,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD { private ulong _divisorLimit; private int _batchSize = 1_024; + private ulong _minK = 1UL; [ThreadStatic] private static GpuUInt128WorkSet _divisorScanGpuWorkSet; @@ -25,6 +26,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD set => _batchSize = Math.Max(1, value); } + public ulong MinK + { + get => _minK; + set => _minK = value < 1UL ? 1UL : value; + } + public void ConfigureFromMaxPrime(ulong maxPrime) { _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); @@ -61,6 +68,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD gpu, prime, allowedMax, + _minK, out processedAll, out divisor); @@ -100,9 +108,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong allowedMax, + ulong minK, out bool processedAll, out ulong foundDivisor) { + foundDivisor = 0UL; + processedAll = true; + // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. // Keeping the guard commented out documents the reasoning for benchmarks and tests. // if (allowedMax < 3UL) @@ -125,6 +137,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD divisor.High = step.High; divisor.Low = step.Low; divisor.Add(1UL); + if (minK > 1UL) + { + GpuUInt128 offset = new(minK - 1UL); + divisor = step * offset; + divisor.Add(1UL); + } + if (divisor.CompareTo(limit) > 0) { processedAll = true; @@ -165,31 +184,66 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD step7 = (byte)(stepLow % 7UL); step11 = (byte)(stepLow % 11UL); - remainder10 = (byte)(divisorLow % 10UL); - remainder8 = (byte)(divisorLow & 7UL); - remainder3 = (byte)(divisorLow % 3UL); - remainder7 = (byte)(divisorLow % 7UL); - remainder11 = (byte)(divisorLow % 11UL); + ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL; + if (maxK == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } - return CheckDivisors64( - gpu, - prime, - stepLow, - limit.Low, - divisorLow, - decimalMask, - step10, - step8, - step3, - step7, - step11, - remainder10, - remainder8, - remainder3, - remainder7, - remainder11, - out processedAll, - out foundDivisor); + ulong startK = minK < 1UL ? 1UL : minK; + bool processedTop = true; + bool processedBottom = true; + bool composite = false; + + if (startK <= maxK) + { + composite = CheckDivisors64Range( + gpu, + prime, + stepLow, + allowedMax, + startK, + maxK, + decimalMask, + step10, + step8, + step3, + step7, + step11, + out processedTop, + out foundDivisor); + if (composite) + { + processedAll = true; + return true; + } + } + + ulong lowerEnd = startK > 1UL ? startK - 1UL : 0UL; + if (lowerEnd >= 1UL) + { + lowerEnd = Math.Min(lowerEnd, maxK); + composite = CheckDivisors64Range( + gpu, + prime, + stepLow, + allowedMax, + 1UL, + lowerEnd, + decimalMask, + step10, + step8, + step3, + step7, + step11, + out processedBottom, + out foundDivisor); + } + + processedAll = processedTop && processedBottom; + return composite; } ulong divisorHigh = divisor.High; @@ -265,6 +319,75 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberD return false; } + private static bool CheckDivisors64Range( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong step, + ulong allowedMax, + ulong startK, + ulong endK, + ushort decimalMask, + byte step10, + byte step8, + byte step3, + byte step7, + byte step11, + out bool processedAll, + out ulong foundDivisor) + { + if (startK < 1UL || endK < startK) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + UInt128 step128 = step; + UInt128 allowedMax128 = allowedMax; + UInt128 startDivisor128 = (step128 * startK) + UInt128.One; + if (startDivisor128 > allowedMax128) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + UInt128 rangeLimit128 = (step128 * endK) + UInt128.One; + if (rangeLimit128 > allowedMax128) + { + rangeLimit128 = allowedMax128; + } + + ulong startDivisor = (ulong)startDivisor128; + ulong rangeLimit = (ulong)rangeLimit128; + + byte remainder10 = (byte)(startDivisor % 10UL); + byte remainder8 = (byte)(startDivisor & 7UL); + byte remainder3 = (byte)(startDivisor % 3UL); + byte remainder7 = (byte)(startDivisor % 7UL); + byte remainder11 = (byte)(startDivisor % 11UL); + + return CheckDivisors64( + gpu, + prime, + step, + rangeLimit, + startDivisor, + decimalMask, + step10, + step8, + step3, + step7, + step11, + remainder10, + remainder8, + remainder3, + remainder7, + remainder11, + out processedAll, + out foundDivisor); + } + private static bool CheckDivisors64( PrimeOrderCalculatorAccelerator gpu, ulong prime, diff --git a/PerfectNumbers.Core/Gpu/AcceleratorPool.cs b/PerfectNumbers.Core/Gpu/AcceleratorPool.cs index 0024155..4cfed88 100644 --- a/PerfectNumbers.Core/Gpu/AcceleratorPool.cs +++ b/PerfectNumbers.Core/Gpu/AcceleratorPool.cs @@ -1,9 +1,10 @@ using ILGPU; using ILGPU.Runtime; +using Open.Collections; namespace PerfectNumbers.Core.Gpu { - public struct AcceleratorPool + public struct AcceleratorPool : IDisposable { internal readonly Context Context; internal readonly Device Device; @@ -28,6 +29,23 @@ namespace PerfectNumbers.Core.Gpu // Intentionally left empty - there's nothing to do here. We're not really renting anything. } + public readonly void Dispose() + { + Accelerator[] accelerators = Accelerators; + HashSet contexts = new(accelerators.Length); + foreach (var accelerator in accelerators) + { + contexts.Add(accelerator.Context); + accelerator.Synchronize(); + accelerator.Dispose(); + } + + foreach (var context in contexts) + { + context.Dispose(); + } + } + public static readonly AcceleratorPool Shared = new(PerfectNumberConstants.RollingAccelerators); private readonly int Capacity; diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index faac3bd..98454a4 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs @@ -1,3 +1,4 @@ +using System.Collections.Concurrent; using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; @@ -9,8 +10,8 @@ public sealed class PrimeOrderCalculatorAccelerator { #region Pool - [ThreadStatic] - private static Queue? _pool; + // [ThreadStatic] + private static ConcurrentQueue? _pool; [MethodImpl(MethodImplOptions.AggressiveInlining)] public static PrimeOrderCalculatorAccelerator Rent(int primeTesterCapacity) @@ -34,7 +35,7 @@ public sealed class PrimeOrderCalculatorAccelerator internal static void Clear() { - Queue? pool = _pool; + var pool = _pool; if (pool != null) { while (pool.TryDequeue(out var gpu)) @@ -46,7 +47,7 @@ public sealed class PrimeOrderCalculatorAccelerator internal static void DisposeAll() { - Queue? pool = _pool; + var pool = _pool; if (pool != null) { while (pool.TryDequeue(out var gpu)) @@ -152,7 +153,7 @@ public sealed class PrimeOrderCalculatorAccelerator private sealed class CalculatorKernels(Accelerator accelerator) { - public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(PrimeOrderKernels.CalculateOrderKernel)).CreateLauncherDelegate, ArrayView1D, int, OrderKernelBuffers>>(); + // public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(PrimeOrderKernels.CalculateOrderKernel)).CreateLauncherDelegate, ArrayView1D, int, OrderKernelBuffers>>(); public readonly Kernel CheckFactorsKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, Stride1D.Dense>, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(PrimeOrderKernels.CheckFactorsKernel)); @@ -162,15 +163,15 @@ public sealed class PrimeOrderCalculatorAccelerator public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> KeepMontgomeryKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery)).CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); - public readonly Action, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> PartialFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(PrimeOrderKernels.PartialFactorKernel)).CreateLauncherDelegate, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + // public readonly Action, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> PartialFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(PrimeOrderKernels.PartialFactorKernel)).CreateLauncherDelegate, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); - public readonly Kernel PollardRhoKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, ArrayView1D, ArrayView1D>(Pow2MontgomeryKernels.TryPollardRhoKernel)); + // public readonly Kernel PollardRhoKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, ArrayView1D, ArrayView1D>(Pow2MontgomeryKernels.TryPollardRhoKernel)); public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> Pow2ModKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(PrimeOrderKernels.Pow2ModKernel)).CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); public readonly Action, GpuUInt128, ArrayView1D> Pow2ModWideKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(PrimeOrderKernels.Pow2ModKernelWide)).CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); - public readonly Action, ArrayView> SharesFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)).CreateLauncherDelegate, ArrayView>>(); + // public readonly Action, ArrayView> SharesFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)).CreateLauncherDelegate, ArrayView>>(); public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan)).CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); @@ -437,16 +438,16 @@ public sealed class PrimeOrderCalculatorAccelerator var kernels = _kernels[acceleratorIndex]; - CalculateOrderKernelLauncher = kernels.CalculateOrderKernelLauncher; + // CalculateOrderKernelLauncher = kernels.CalculateOrderKernelLauncher; CheckFactorsKernel = kernels.CheckFactorsKernel; ConvertToStandardKernelLauncher = kernels.ConvertToStandardKernelLauncher; HeuristicCombinedTrialDivisionKernelLauncher = kernels.HeuristicCombinedTrialDivisionKernelLauncher; KeepMontgomeryKernelLauncher = kernels.KeepMontgomeryKernelLauncher; - PartialFactorKernelLauncher = kernels.PartialFactorKernelLauncher; - PollardRhoKernel = kernels.PollardRhoKernel; + // PartialFactorKernelLauncher = kernels.PartialFactorKernelLauncher; + // PollardRhoKernel = kernels.PollardRhoKernel; Pow2ModKernelLauncher = kernels.Pow2ModKernelLauncher; Pow2ModWideKernelLauncher = kernels.Pow2ModWideKernelLauncher; - SharesFactorKernelLauncher = kernels.SharesFactorKernelLauncher; + // SharesFactorKernelLauncher = kernels.SharesFactorKernelLauncher; SmallPrimeFactorKernelLauncher = kernels.SmallPrimeFactorKernelLauncher; SmallPrimeSieveKernelLauncher = kernels.SmallPrimeSieveKernelLauncher; SpecialMaxKernelLauncher = kernels.SpecialMaxKernelLauncher; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 5b2590c..b28e598 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -18,6 +18,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN // so the synchronization fields from the previous implementation remain commented out here. // private readonly object _sync = new(); private ulong _divisorLimit; + private ulong _minK = 1UL; // private bool _isConfigured; private readonly ConcurrentBag[] _resourcePool = [.. AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentBag())]; @@ -36,6 +37,12 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN set => GpuBatchSize = value; } + public ulong MinK + { + get => _minK; + set => _minK = value < 1UL ? 1UL : value; + } + public void ConfigureFromMaxPrime(ulong maxPrime) { // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, @@ -95,6 +102,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN gpu, prime, allowedMax, + _minK, resources.CheckDivisorKernel, resources.DivisorDataBuffer, resources.OffsetBuffer, @@ -162,6 +170,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong allowedMax, + ulong minK, Kernel kernel, MemoryBuffer1D divisorDataBuffer, MemoryBuffer1D offsetBuffer, @@ -183,8 +192,8 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN int batchCapacity = (int)divisorDataBuffer.Length; bool composite = false; bool processedAll; - processedCount = 0UL; - lastProcessed = 0UL; + ulong processedCountLocal = 0UL; + ulong lastProcessedLocal = 0UL; var acceleratorIndex = gpu.AcceleratorIndex; // EvenPerfectBitScanner always configures non-empty GPU divisor buffers, so the defensive guard stays @@ -208,11 +217,13 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN if (!hasCandidates) { coveredRange = true; + processedCount = processedCountLocal; + lastProcessed = lastProcessedLocal; return false; } ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; - ulong remainingCount = maxK; + ulong startK = minK < 1UL ? 1UL : minK; ulong stepValue = (ulong)twoP128; byte step10 = (byte)(stepValue % 10UL); @@ -221,13 +232,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN byte step7 = (byte)(stepValue % 7UL); byte step11 = (byte)(stepValue % 11UL); - UInt128 currentDivisor128 = firstDivisor128; - ulong currentDivisorValue = (ulong)currentDivisor128; - byte remainder10 = (byte)(currentDivisorValue % 10UL); - byte remainder8 = (byte)(currentDivisorValue & 7UL); - byte remainder3 = (byte)(currentDivisorValue % 3UL); - byte remainder7 = (byte)(currentDivisorValue % 7UL); - byte remainder11 = (byte)(currentDivisorValue % 11UL); LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); @@ -243,11 +247,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN Span offsetSpan = offsetStorage; Span countSpan = countStorage; Span cycleSpan = cycleStorage; - ref GpuDivisorPartialData divisorDataRef = ref MemoryMarshal.GetReference(divisorDataSpan); - ref int offsetRef = ref MemoryMarshal.GetReference(offsetSpan); - ref int countRef = ref MemoryMarshal.GetReference(countSpan); - ref ulong cycleRef = ref MemoryMarshal.GetReference(cycleSpan); - ref ulong exponentRef = ref MemoryMarshal.GetReference(exponentSpan); ArrayView1D divisorDataView = divisorDataBuffer.View; ArrayView1D offsetView = offsetBuffer.View; @@ -257,100 +256,258 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN ArrayView1D hitsView = hitsBuffer.View; ArrayView1D hitIndexView = hitIndexBuffer.View; Queue divisorPool = MontgomeryDivisorDataPool.Shared; - while (remainingCount > 0UL && !composite) + static bool ProcessRange( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + UInt128 twoP128, + UInt128 allowedMax128, + byte step10, + byte step8, + byte step3, + byte step7, + byte step11, + ushort decimalMask, + ArrayView1D divisorDataView, + ArrayView1D offsetView, + ArrayView1D countView, + ArrayView1D cycleView, + ArrayView1D exponentViewDevice, + ArrayView1D hitsView, + ArrayView1D hitIndexView, + Span divisorSpan, + Span exponentSpan, + Span divisorDataSpan, + Span offsetSpan, + Span countSpan, + Span cycleSpan, + Queue divisorPool, + int chunkCountBaseline, + int acceleratorIndex, + ref bool composite, + ref ulong processedCountLocal, + ref ulong lastProcessedLocal, + ulong rangeStartK, + ulong rangeEndK, + Kernel kernel) { - int chunkCount = chunkCountBaseline; - if ((ulong)chunkCount > remainingCount) + if (rangeStartK < 1UL || rangeEndK < rangeStartK) { - chunkCount = (int)remainingCount; + return true; } - UInt128 nextDivisor128 = currentDivisor128; - int admissibleCount = 0; - - byte localRemainder10 = remainder10; - byte localRemainder8 = remainder8; - byte localRemainder3 = remainder3; - byte localRemainder7 = remainder7; - byte localRemainder11 = remainder11; - - for (int i = 0; i < chunkCount; i++) + UInt128 startDivisor128 = (twoP128 * rangeStartK) + UInt128.One; + if (startDivisor128 > allowedMax128) { - bool passesSmallModuli = localRemainder3 != 0 && localRemainder7 != 0 && localRemainder11 != 0; - if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0) - { - ulong divisorValue = (ulong)nextDivisor128; - MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue); - ulong divisorCycle = ResolveDivisorCycle(gpu, divisorValue, prime, in montgomeryData); - if (divisorCycle == prime) - { - processedCount += (ulong)(i + 1); - lastProcessed = divisorValue; - coveredRange = true; - return true; - } + return true; + } - divisorPool.Return(montgomeryData); - divisorSpan[admissibleCount] = divisorValue; - divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); - offsetSpan[admissibleCount] = admissibleCount; - countSpan[admissibleCount] = 1; - cycleSpan[admissibleCount] = divisorCycle; - exponentSpan[admissibleCount] = prime; - admissibleCount++; + UInt128 rangeLimit128 = (twoP128 * rangeEndK) + UInt128.One; + if (rangeLimit128 > allowedMax128) + { + rangeLimit128 = allowedMax128; + } + + UInt128 rangeCount128 = ((rangeLimit128 - startDivisor128) / twoP128) + UInt128.One; + ulong rangeRemaining = rangeCount128 > ulong.MaxValue ? ulong.MaxValue : (ulong)rangeCount128; + UInt128 currentRangeDivisor128 = startDivisor128; + ulong currentRangeDivisorValue = (ulong)currentRangeDivisor128; + byte rangeRemainder10 = (byte)(currentRangeDivisorValue % 10UL); + byte rangeRemainder8 = (byte)(currentRangeDivisorValue & 7UL); + byte rangeRemainder3 = (byte)(currentRangeDivisorValue % 3UL); + byte rangeRemainder7 = (byte)(currentRangeDivisorValue % 7UL); + byte rangeRemainder11 = (byte)(currentRangeDivisorValue % 11UL); + + while (rangeRemaining > 0UL && !composite) + { + int chunkCount = chunkCountBaseline; + if ((ulong)chunkCount > rangeRemaining) + { + chunkCount = (int)rangeRemaining; } - nextDivisor128 += twoP128; - localRemainder10 = AddMod10(localRemainder10, step10); - localRemainder8 = AddMod8(localRemainder8, step8); - localRemainder3 = AddMod3(localRemainder3, step3); - localRemainder7 = AddMod7(localRemainder7, step7); - localRemainder11 = AddMod11(localRemainder11, step11); + UInt128 nextDivisor128 = currentRangeDivisor128; + int admissibleCount = 0; + + byte localRemainder10 = rangeRemainder10; + byte localRemainder8 = rangeRemainder8; + byte localRemainder3 = rangeRemainder3; + byte localRemainder7 = rangeRemainder7; + byte localRemainder11 = rangeRemainder11; + + for (int i = 0; i < chunkCount; i++) + { + bool passesSmallModuli = localRemainder3 != 0 && localRemainder7 != 0 && localRemainder11 != 0; + if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0) + { + ulong divisorValue = (ulong)nextDivisor128; + MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue); + ulong divisorCycle = ResolveDivisorCycle(gpu, divisorValue, prime, in montgomeryData); + if (divisorCycle == prime) + { + processedCountLocal += (ulong)(i + 1); + lastProcessedLocal = divisorValue; + composite = true; + return false; + } + + divisorPool.Return(montgomeryData); + divisorSpan[admissibleCount] = divisorValue; + divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); + offsetSpan[admissibleCount] = admissibleCount; + countSpan[admissibleCount] = 1; + cycleSpan[admissibleCount] = divisorCycle; + exponentSpan[admissibleCount] = prime; + admissibleCount++; + } + + nextDivisor128 += twoP128; + localRemainder10 = AddMod10(localRemainder10, step10); + localRemainder8 = AddMod8(localRemainder8, step8); + localRemainder3 = AddMod3(localRemainder3, step3); + localRemainder7 = AddMod7(localRemainder7, step7); + localRemainder11 = AddMod11(localRemainder11, step11); + } + + rangeRemainder10 = localRemainder10; + rangeRemainder8 = localRemainder8; + rangeRemainder3 = localRemainder3; + rangeRemainder7 = localRemainder7; + rangeRemainder11 = localRemainder11; + + processedCountLocal += (ulong)chunkCount; + + UInt128 lastDivisor128 = nextDivisor128 - twoP128; + lastProcessedLocal = (ulong)lastDivisor128; + currentRangeDivisor128 = nextDivisor128; + rangeRemaining -= (ulong)chunkCount; + + if (admissibleCount == 0) + { + continue; + } + + ref GpuDivisorPartialData divisorDataRef = ref MemoryMarshal.GetReference(divisorDataSpan); + ref int offsetRef = ref MemoryMarshal.GetReference(offsetSpan); + ref int countRef = ref MemoryMarshal.GetReference(countSpan); + ref ulong cycleRef = ref MemoryMarshal.GetReference(cycleSpan); + ref ulong exponentRef = ref MemoryMarshal.GetReference(exponentSpan); + + var stream = AcceleratorStreamPool.Rent(acceleratorIndex); + divisorDataView.CopyFromCPU(stream, ref divisorDataRef, admissibleCount); + offsetView.CopyFromCPU(stream, ref offsetRef, admissibleCount); + countView.CopyFromCPU(stream, ref countRef, admissibleCount); + cycleView.CopyFromCPU(stream, ref cycleRef, admissibleCount); + exponentViewDevice.CopyFromCPU(stream, ref exponentRef, admissibleCount); + + int sentinel = int.MaxValue; + hitIndexView.CopyFromCPU(stream, ref sentinel, 1); + + kernel.Launch(stream, 1, new Index1D(admissibleCount), divisorDataView, offsetView, countView, exponentViewDevice, cycleView, hitsView, hitIndexView); + + int firstHit = sentinel; + hitIndexView.CopyToCPU(stream, ref firstHit, 1); + stream.Synchronize(); + AcceleratorStreamPool.Return(acceleratorIndex, stream); + int hitIndex = firstHit >= admissibleCount ? -1 : firstHit; + bool hitFound = hitIndex >= 0; + composite = hitFound; + int lastIndex = admissibleCount - 1; + lastProcessedLocal = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; } - remainder10 = localRemainder10; - remainder8 = localRemainder8; - remainder3 = localRemainder3; - remainder7 = localRemainder7; - remainder11 = localRemainder11; - - processedCount += (ulong)chunkCount; - - UInt128 lastDivisor128 = nextDivisor128 - twoP128; - lastProcessed = (ulong)lastDivisor128; - currentDivisor128 = nextDivisor128; - remainingCount -= (ulong)chunkCount; - - if (admissibleCount == 0) - { - continue; - } - - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - divisorDataView.CopyFromCPU(stream, ref divisorDataRef, admissibleCount); - offsetView.CopyFromCPU(stream, ref offsetRef, admissibleCount); - countView.CopyFromCPU(stream, ref countRef, admissibleCount); - cycleView.CopyFromCPU(stream, ref cycleRef, admissibleCount); - exponentViewDevice.CopyFromCPU(stream, ref exponentRef, admissibleCount); - - int sentinel = int.MaxValue; - hitIndexView.CopyFromCPU(stream, ref sentinel, 1); - - kernel.Launch(stream, 1, new Index1D(admissibleCount), divisorDataView, offsetView, countView, exponentViewDevice, cycleView, hitsView, hitIndexView); - - int firstHit = sentinel; - hitIndexView.CopyToCPU(stream, ref firstHit, 1); - stream.Synchronize(); - AcceleratorStreamPool.Return(acceleratorIndex, stream); - int hitIndex = firstHit >= admissibleCount ? -1 : firstHit; - bool hitFound = hitIndex >= 0; - composite = hitFound; - int lastIndex = admissibleCount - 1; - lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; + return rangeRemaining == 0UL; } - processedAll = remainingCount == 0UL; - coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); + bool processedTop = true; + bool processedBottom = true; + + if (startK <= maxK) + { + processedTop = ProcessRange( + gpu, + prime, + twoP128, + allowedMax128, + step10, + step8, + step3, + step7, + step11, + decimalMask, + divisorDataView, + offsetView, + countView, + cycleView, + exponentViewDevice, + hitsView, + hitIndexView, + divisorSpan, + exponentSpan, + divisorDataSpan, + offsetSpan, + countSpan, + cycleSpan, + divisorPool, + chunkCountBaseline, + acceleratorIndex, + ref composite, + ref processedCountLocal, + ref lastProcessedLocal, + startK, + maxK, + kernel); + if (composite) + { + coveredRange = true; + processedCount = processedCountLocal; + lastProcessed = lastProcessedLocal; + return true; + } + } + + ulong lowerEnd = startK > 1UL ? startK - 1UL : 0UL; + if (!composite && lowerEnd >= 1UL) + { + lowerEnd = Math.Min(lowerEnd, maxK); + processedBottom = ProcessRange( + gpu, + prime, + twoP128, + allowedMax128, + step10, + step8, + step3, + step7, + step11, + decimalMask, + divisorDataView, + offsetView, + countView, + cycleView, + exponentViewDevice, + hitsView, + hitIndexView, + divisorSpan, + exponentSpan, + divisorDataSpan, + offsetSpan, + countSpan, + cycleSpan, + divisorPool, + chunkCountBaseline, + acceleratorIndex, + ref composite, + ref processedCountLocal, + ref lastProcessedLocal, + 1UL, + lowerEnd, + kernel); + } + + processedAll = processedTop && processedBottom; + coveredRange = composite || processedAll; + processedCount = processedCountLocal; + lastProcessed = lastProcessedLocal; return composite; } diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index b060f43..e8edb89 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -6,6 +6,8 @@ public interface IMersenneNumberDivisorByDivisorTester { int BatchSize { get; set; } + ulong MinK { get; set; } + void ConfigureFromMaxPrime(ulong maxPrime); ulong DivisorLimit { get; } diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 2033992..903ed1c 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -8,7 +8,7 @@ public static class PerfectNumberConstants public const int DefaultSmallPrimeFactorSlotCount = 64; public const int DefaultSpecialMaxFactorCapacity = 1024; public const ulong ExtraDivisorCycleSearchLimit = 64UL; - public const int GpuFrequency = 4 - 1; + public const int GpuFrequency = 256 - 1; public const int MaxQForDivisorCycles = 4_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxOddPowersCount = 128; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 42e1334..1e59cc4 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -54,11 +54,11 @@ internal static partial class PrimeOrderCalculator ulong phi = prime - 1UL; - if (PrimeOrderGpuHeuristics.TryCalculateOrder(gpu, prime, previousOrder, config, divisorData, out ulong gpuOrder)) + // if (PrimeOrderGpuHeuristics.TryCalculateOrder(gpu, prime, previousOrder, config, divisorData, out ulong gpuOrder)) // if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(gpu, prime, previousOrder, config, divisorData, out ulong gpuOrder)) - { - return gpuOrder; - } + // { + // return gpuOrder; + // } PartialFactorResult phiFactors = PartialFactor(gpu, phi, config); @@ -137,7 +137,7 @@ internal static partial class PrimeOrderCalculator Span stackBuffer = stackalloc ulong[length]; - if (length <= 8) + if (length <= 7) { return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData); } @@ -1106,19 +1106,19 @@ internal static partial class PrimeOrderCalculator uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; ulong remaining = value; - bool gpuFactored = false; - if (IsGpuHeuristicDevice) - { - gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out factorCount, out remaining); - } + // bool gpuFactored = false; + // if (IsGpuHeuristicDevice) + // { + // gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out factorCount, out remaining); + // } List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); Stack? compositeStack = null; PrimeTester primeTester = PrimeTester.Exclusive; // HeuristicPrimeTester? tester = _tester ??= new(); - if (!gpuFactored) - { + // if (!gpuFactored) + // { PopulateSmallPrimeFactorsCpu( value, limit, @@ -1128,7 +1128,7 @@ internal static partial class PrimeOrderCalculator ref exponentSlotArray, out factorCount, out remaining); - } + // } if (remaining > 1UL) @@ -1176,7 +1176,7 @@ internal static partial class PrimeOrderCalculator // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeCpu(composite); // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); bool runOnGpu = RunOnGpu(); - bool isPrime = !runOnGpu ? HeuristicCombinedPrimeTester.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + bool isPrime = !runOnGpu ? PrimeTester.IsPrime(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); // bool isPrime = PrimeTester.IsPrime(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); @@ -1244,15 +1244,15 @@ internal static partial class PrimeOrderCalculator { // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeCpu(composite); // bool isPrime = PrimeTester.IsPrime(composite); - // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(composite); + bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); // bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); // Atomic.Add(ref _partialFactorPendingHits, 1UL); // Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}"); - bool isPrime = RunOnGpu() - ? HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite) - : HeuristicCombinedPrimeTester.IsPrimeCpu(composite); + // bool isPrime = RunOnGpu() + // ? HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite) + // : PrimeTester.IsPrime(composite); // bool isPrime = PrimeTester.IsPrime(composite); // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); @@ -1287,7 +1287,7 @@ internal static partial class PrimeOrderCalculator { // cofactorIsPrime = PrimeTester.IsPrime(cofactor); // cofactorIsPrime = HeuristicCombinedPrimeTester.IsPrimeCpu(cofactor); - // cofactorIsPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(cofactor); + cofactorIsPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor); // cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None); // Atomic.Add(ref _partialFactorCofactorHits, 1UL); @@ -1295,9 +1295,9 @@ internal static partial class PrimeOrderCalculator // cofactorIsPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor); // cofactorIsPrime = PrimeTester.IsPrimeGpu(gpu, cofactor); - cofactorIsPrime = RunOnGpu() - ? HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor) - : HeuristicCombinedPrimeTester.IsPrimeCpu(cofactor); + // cofactorIsPrime = RunOnGpu() + // ? HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor) + // : PrimeTester.IsPrime(cofactor); // cofactorIsPrime = PrimeTester.IsPrime(cofactor); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1385,7 +1385,7 @@ internal static partial class PrimeOrderCalculator private static bool RunOnGpu() { // return false; - return true; + // return true; // int cpuCount; // lock(_sync)