From 94a4ed83283a5386b149b6ff04eff91756433129 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 27 Oct 2025 21:24:43 +0100 Subject: [PATCH 1/2] Hoist prime order thresholds in CPU heuristics --- EvenPerfectBitScanner/CliArguments.cs | 79 +++-- .../IO/CalculationResultsFile.cs | 69 +++- ...MersenneNumberDivisorByDivisorCpuTester.cs | 334 ++++++++++++++---- PerfectNumbers.Core/DivisorGenerator.cs | 28 +- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 53 ++- ...MersenneNumberDivisorByDivisorGpuTester.cs | 273 ++++++++------ PerfectNumbers.Core/MersenneDivisorCycles.cs | 60 ++-- .../MersenneNumberDivisorByDivisorTester.cs | 10 +- .../PrimeOrderCalculator.Cpu.cs | 77 ++-- docs/ByDivisorPathMap.md | 35 ++ docs/HeuristicIsPrimePlan.md | 2 +- todo_list.txt | 1 - 12 files changed, 670 insertions(+), 351 deletions(-) create mode 100644 docs/ByDivisorPathMap.md diff --git a/EvenPerfectBitScanner/CliArguments.cs b/EvenPerfectBitScanner/CliArguments.cs index 7c9b240..3145ac4 100644 --- a/EvenPerfectBitScanner/CliArguments.cs +++ b/EvenPerfectBitScanner/CliArguments.cs @@ -202,14 +202,15 @@ internal readonly struct CliArguments if (argument.StartsWith("--prime=", StringComparison.OrdinalIgnoreCase)) { - startPrime = Utf8CliParser.ParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1)); + ReadOnlySpan primeValue = argument.AsSpan("--prime=".Length); + startPrime = Utf8CliParser.ParseUInt64(primeValue); startPrimeProvided = true; continue; } if (argument.StartsWith("--max-prime=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan(argument.IndexOf('=') + 1); + ReadOnlySpan value = argument.AsSpan("--max-prime=".Length); if (UInt128.TryParse(value, NumberStyles.None, CultureInfo.InvariantCulture, out UInt128 parsedMaxPrime)) { maxPrimeLimit = parsedMaxPrime; @@ -238,13 +239,13 @@ internal readonly struct CliArguments if (argument.StartsWith("--threads=", StringComparison.OrdinalIgnoreCase)) { - threadCount = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + threadCount = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--threads=".Length))); continue; } if (argument.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan(argument.IndexOf('=') + 1); + ReadOnlySpan value = argument.AsSpan("--mersenne=".Length); if (value.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) { kernelType = GpuKernelType.Pow2Mod; @@ -281,7 +282,7 @@ internal readonly struct CliArguments if (argument.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedLimit)) + if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--divisor-cycles-limit=".Length), out ulong parsedLimit)) { divisorCyclesSearchLimit = parsedLimit; } @@ -291,7 +292,7 @@ internal readonly struct CliArguments if (argument.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedResidueMax)) + if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--residue-max-k=".Length), out ulong parsedResidueMax)) { residueKMax = parsedResidueMax; } @@ -301,7 +302,7 @@ internal readonly struct CliArguments if (argument.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase)) { - useMersenneOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + useMersenneOnGpu = !argument.AsSpan("--mersenne-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } @@ -319,7 +320,7 @@ internal readonly struct CliArguments if (argument.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedLimit)) + if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--order-warmup-limit=".Length), out ulong parsedLimit)) { orderWarmupLimitOverride = parsedLimit; } @@ -335,14 +336,14 @@ internal readonly struct CliArguments if (argument.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan(argument.IndexOf('=') + 1); + ReadOnlySpan value = argument.AsSpan("--ntt=".Length); nttBackend = value.Equals("staged", StringComparison.OrdinalIgnoreCase) ? NttBackend.Staged : NttBackend.Reference; continue; } if (argument.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan(argument.IndexOf('=') + 1); + ReadOnlySpan value = argument.AsSpan("--mod-reduction=".Length); if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase)) { modReductionMode = ModReductionMode.Mont64; @@ -365,25 +366,25 @@ internal readonly struct CliArguments if (argument.StartsWith("--primes-device=", StringComparison.OrdinalIgnoreCase)) { - forcePrimeKernelsOnCpu = argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + forcePrimeKernelsOnCpu = argument.AsSpan("--primes-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } if (argument.StartsWith("--order-device=", StringComparison.OrdinalIgnoreCase)) { - useOrderOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + useOrderOnGpu = !argument.AsSpan("--order-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } if (argument.StartsWith("--rle-blacklist=", StringComparison.OrdinalIgnoreCase)) { - rleBlacklistPath = argument[(argument.IndexOf('=') + 1)..]; + rleBlacklistPath = argument["--rle-blacklist=".Length..]; continue; } if (argument.StartsWith("--rle-hard-max=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedMax)) + if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--rle-hard-max=".Length), out ulong parsedMax)) { rleHardMaxP = parsedMax; } @@ -393,14 +394,14 @@ internal readonly struct CliArguments if (argument.StartsWith("--rle-only-last7=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan(argument.IndexOf('=') + 1); + ReadOnlySpan value = argument.AsSpan("--rle-only-last7=".Length); rleOnlyLast7 = !value.Equals("false", StringComparison.OrdinalIgnoreCase) && !value.Equals("0", StringComparison.OrdinalIgnoreCase); continue; } if (argument.StartsWith("--zero-hard=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseDouble(argument.AsSpan(argument.IndexOf('=') + 1), out double zeroFraction)) + if (Utf8CliParser.TryParseDouble(argument.AsSpan("--zero-hard=".Length), out double zeroFraction)) { zeroFractionHard = zeroFraction; } @@ -410,7 +411,7 @@ internal readonly struct CliArguments if (argument.StartsWith("--zero-conj=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan(argument.IndexOf('=') + 1); + ReadOnlySpan value = argument.AsSpan("--zero-conj=".Length); int colon = value.IndexOf(':'); if (colon > 0) { @@ -427,87 +428,97 @@ internal readonly struct CliArguments if (argument.StartsWith("--lucas=", StringComparison.OrdinalIgnoreCase)) { - useMersenneOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + useMersenneOnGpu = !argument.AsSpan("--lucas=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } - if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase) || - argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase) || - argument.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase)) + if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase)) { - forcePrimeKernelsOnCpu = argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + forcePrimeKernelsOnCpu = argument.AsSpan("--primes=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); + continue; + } + + if (argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase)) + { + forcePrimeKernelsOnCpu = argument.AsSpan("--gpu-kernels=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); + continue; + } + + if (argument.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase)) + { + forcePrimeKernelsOnCpu = argument.AsSpan("--accelerator=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } if (argument.StartsWith("--results-dir=", StringComparison.OrdinalIgnoreCase)) { - resultsDirectory = argument[(argument.IndexOf('=') + 1)..]; + resultsDirectory = argument["--results-dir=".Length..]; continue; } if (argument.StartsWith("--results-prefix=", StringComparison.OrdinalIgnoreCase)) { - resultsPrefix = argument[(argument.IndexOf('=') + 1)..]; + resultsPrefix = argument["--results-prefix=".Length..]; continue; } if (argument.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase)) { - gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-threads=".Length))); continue; } if (argument.StartsWith("--gpu-prime-batch=", StringComparison.OrdinalIgnoreCase)) { - gpuPrimeBatch = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + gpuPrimeBatch = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-batch=".Length))); continue; } if (argument.StartsWith("--ll-slice=", StringComparison.OrdinalIgnoreCase)) { - sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--ll-slice=".Length))); continue; } if (argument.StartsWith("--gpu-scan-batch=", StringComparison.OrdinalIgnoreCase)) { - scanBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + scanBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-scan-batch=".Length))); continue; } if (argument.StartsWith("--block-size=", StringComparison.OrdinalIgnoreCase)) { - blockSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + blockSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--block-size=".Length))); continue; } if (argument.StartsWith("--filter-p=", StringComparison.OrdinalIgnoreCase)) { - filterFile = argument[(argument.IndexOf('=') + 1)..]; + filterFile = argument["--filter-p=".Length..]; continue; } if (argument.StartsWith("--write-batch-size=", StringComparison.OrdinalIgnoreCase)) { - writeBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + writeBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--write-batch-size=".Length))); continue; } if (argument.StartsWith("--divisor-cycles=", StringComparison.OrdinalIgnoreCase)) { - cyclesPath = argument[(argument.IndexOf('=') + 1)..]; + cyclesPath = argument["--divisor-cycles=".Length..]; continue; } if (argument.StartsWith("--divisor-cycles-device=", StringComparison.OrdinalIgnoreCase)) { - useGpuCycles = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + useGpuCycles = !argument.AsSpan("--divisor-cycles-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } if (argument.StartsWith("--divisor-cycles-batch=", StringComparison.OrdinalIgnoreCase)) { - cyclesBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1))); + cyclesBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--divisor-cycles-batch=".Length))); continue; } diff --git a/EvenPerfectBitScanner/IO/CalculationResultsFile.cs b/EvenPerfectBitScanner/IO/CalculationResultsFile.cs index e716002..18fbbdb 100644 --- a/EvenPerfectBitScanner/IO/CalculationResultsFile.cs +++ b/EvenPerfectBitScanner/IO/CalculationResultsFile.cs @@ -14,7 +14,7 @@ internal static class CalculationResultsFile using FileStream readStream = new(candidateFile, FileMode.Open, FileAccess.Read, FileShare.ReadWrite); using StreamReader reader = new(readStream); string? line; - ReadOnlySpan span; + ReadOnlySpan trimmed; while ((line = reader.ReadLine()) is not null) { if (string.IsNullOrWhiteSpace(line)) @@ -22,36 +22,67 @@ internal static class CalculationResultsFile continue; } - span = line.AsSpan().Trim(); - if (span.IsEmpty || span[0] == '#') + trimmed = line.AsSpan().Trim(); + if (trimmed.IsEmpty || trimmed[0] == '#') { continue; } + + int length = trimmed.Length; int index = 0; - while (index < span.Length) + if (!maxPrimeConfigured) { - if (!char.IsDigit(span[index])) + while (index < length) { - index++; - continue; - } + char current = trimmed[index]; + uint digit = (uint)(current - '0'); + if (digit > 9U) + { + index++; + continue; + } - int start = index; - index++; - while (index < span.Length && char.IsDigit(span[index])) - { - index++; - } + int start = index++; + while (index < length && (uint)(trimmed[index] - '0') <= 9U) + { + index++; + } - if (Utf8CliParser.TryParseUInt64(span[start..index], out ulong parsed)) - { - if (!maxPrimeConfigured || (UInt128)parsed <= maxPrimeLimit) + if (Utf8CliParser.TryParseUInt64(trimmed[start..index], out ulong parsed)) { candidates.Add(parsed); } - else + } + } + else + { + UInt128 limit = maxPrimeLimit; + while (index < length) + { + char current = trimmed[index]; + uint digit = (uint)(current - '0'); + if (digit > 9U) { - skippedByLimit++; + index++; + continue; + } + + int start = index++; + while (index < length && (uint)(trimmed[index] - '0') <= 9U) + { + index++; + } + + if (Utf8CliParser.TryParseUInt64(trimmed[start..index], out ulong parsed)) + { + if ((UInt128)parsed <= limit) + { + candidates.Add(parsed); + } + else + { + skippedByLimit++; + } } } } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 5254d7a..e68668d 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -7,10 +7,7 @@ namespace PerfectNumbers.Core.Cpu; public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester { - private readonly object _sync = new(); - private ulong _divisorLimit; - private ulong _lastStatusDivisor; private int _batchSize = 1_024; [ThreadStatic] @@ -32,7 +29,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv public void ConfigureFromMaxPrime(ulong maxPrime) { _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - _lastStatusDivisor = 0UL; } public ulong DivisorLimit @@ -60,22 +56,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // return true; // } - ulong processedCount; bool processedAll; bool composite = CheckDivisors( prime, allowedMax, - out processedAll, - out processedCount); - - if (processedCount > 0UL) - { - lock (_sync) - { - UpdateStatusUnsafe(processedCount); - } - } + out processedAll); if (composite) { @@ -111,10 +97,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool CheckDivisors( ulong prime, ulong allowedMax, - out bool processedAll, - out ulong processedCount) + out bool processedAll) { - processedCount = 0UL; processedAll = false; // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. @@ -151,44 +135,86 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong stepHigh = step.High; ulong stepLow = step.Low; - byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL); - byte step8 = (byte)(stepLow % 8UL); - byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL); - byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL); - byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL); - byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL); - - ulong divisorHigh = divisor.High; + ulong limitHigh = limit.High; ulong divisorLow = divisor.Low; - byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL); - byte remainder8 = (byte)(divisorLow % 8UL); - byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL); - byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL); - byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL); - byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL); // Keep the divisibility filters aligned with the divisor-cycle generator so the // CPU path never requests cycles that were skipped during cache creation. LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - bool lastIsSeven = lastDigit == LastDigit.Seven; - const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); - const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); + + byte step10; + byte step8; + byte step5; + byte step3; + byte step7; + byte step11; + + byte remainder10; + byte remainder8; + byte remainder5; + byte remainder3; + byte remainder7; + byte remainder11; + + if (stepHigh == 0UL && limitHigh == 0UL) + { + step10 = (byte)(stepLow % 10UL); + step8 = (byte)(stepLow & 7UL); + step5 = (byte)(stepLow % 5UL); + step3 = (byte)(stepLow % 3UL); + step7 = (byte)(stepLow % 7UL); + step11 = (byte)(stepLow % 11UL); + + remainder10 = (byte)(divisorLow % 10UL); + remainder8 = (byte)(divisorLow & 7UL); + remainder5 = (byte)(divisorLow % 5UL); + remainder3 = (byte)(divisorLow % 3UL); + remainder7 = (byte)(divisorLow % 7UL); + remainder11 = (byte)(divisorLow % 11UL); + + return CheckDivisors64( + prime, + stepLow, + limit.Low, + divisorLow, + decimalMask, + step10, + step8, + step5, + step3, + step7, + step11, + remainder10, + remainder8, + remainder5, + remainder3, + remainder7, + remainder11, + out processedAll); + } + + ulong divisorHigh = divisor.High; + step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL); + step8 = (byte)(stepLow % 8UL); + step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL); + step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL); + step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL); + step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL); + + remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL); + remainder8 = (byte)(divisorLow % 8UL); + remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL); + remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL); + remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL); + remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL); while (divisor.CompareTo(limit) <= 0) { - ulong candidate = divisor.Low; - processedCount++; - - ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; - bool admissible = ((decimalMask >> remainder10) & 1) != 0; - - if (admissible - && (remainder8 == 1 || remainder8 == 7) - && remainder3 != 0 - && remainder5 != 0 - && remainder7 != 0 - && remainder11 != 0) + bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0; + if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) { + ulong candidate = divisor.Low; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); ulong divisorCycle; // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. @@ -233,16 +259,141 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } divisor.Add(step); - remainder10 = AddMod(remainder10, step10, (byte)10); - remainder8 = AddMod(remainder8, step8, (byte)8); - remainder5 = AddMod(remainder5, step5, (byte)5); - remainder3 = AddMod(remainder3, step3, (byte)3); - remainder7 = AddMod(remainder7, step7, (byte)7); - remainder11 = AddMod(remainder11, step11, (byte)11); + remainder10 = AddMod10(remainder10, step10); + remainder8 = AddMod8(remainder8, step8); + remainder5 = AddMod5(remainder5, step5); + remainder3 = AddMod3(remainder3, step3); + remainder7 = AddMod7(remainder7, step7); + remainder11 = AddMod11(remainder11, step11); + } + processedAll = true; + return false; + } + + private static bool CheckDivisors64( + ulong prime, + ulong step, + ulong limit, + ulong divisor, + ushort decimalMask, + byte step10, + byte step8, + byte step5, + byte step3, + byte step7, + byte step11, + byte remainder10, + byte remainder8, + byte remainder5, + byte remainder3, + byte remainder7, + byte remainder11, + out bool processedAll) + { + processedAll = false; + + bool canAdvance = step <= limit; + if (!canAdvance) + { + if (divisor <= limit) + { + if (remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + divisor, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneDivisorCycles.CalculateCycleLength( + divisor, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + if (divisorCycle == prime) + { + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + } + + processedAll = true; + return false; } - processedAll = divisor.CompareTo(limit) > 0; - return false; + if (divisor > limit) + { + processedAll = true; + return false; + } + + ulong remainingIterations = ((limit - divisor) / step) + 1UL; + + while (true) + { + bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0; + if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + divisor, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneDivisorCycles.CalculateCycleLength( + divisor, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + if (divisorCycle == prime) + { + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + remainingIterations--; + if (remainingIterations == 0UL) + { + processedAll = true; + return false; + } + + divisor += step; + remainder10 = AddMod10(remainder10, step10); + remainder8 = AddMod8(remainder8, step8); + remainder5 = AddMod5(remainder5, step5); + remainder3 = AddMod3(remainder3, step3); + remainder7 = AddMod7(remainder7, step7); + remainder11 = AddMod11(remainder11, step11); + } } private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) @@ -251,35 +402,80 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return residue == 1UL ? (byte)1 : (byte)0; } - private static byte AddMod(byte value, byte delta, byte modulus) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod3(byte value, byte delta) { + const int Modulus = 3; int sum = value + delta; - if (sum >= modulus) + + if (sum >= Modulus) { - sum -= modulus; + sum -= Modulus; } return (byte)sum; } - private void UpdateStatusUnsafe(ulong processedCount) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod5(byte value, byte delta) { - if (processedCount == 0UL) + const int Modulus = 5; + int sum = value + delta; + + if (sum >= Modulus) { - return; + sum -= Modulus; } - ulong interval = PerfectNumberConstants.ConsoleInterval; - if (interval == 0UL) + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod7(byte value, byte delta) + { + const int Modulus = 7; + int sum = value + delta; + + if (sum >= Modulus) { - _lastStatusDivisor = 0UL; - return; + sum -= Modulus; } - ulong total = _lastStatusDivisor + processedCount; - // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI - // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly. - _lastStatusDivisor = total % interval; + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod10(byte value, byte delta) + { + const int Modulus = 10; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod11(byte value, byte delta) + { + const int Modulus = 11; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod8(byte value, byte delta) + { + return (byte)((value + delta) & 7); } private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) diff --git a/PerfectNumbers.Core/DivisorGenerator.cs b/PerfectNumbers.Core/DivisorGenerator.cs index 0001267..fdb648c 100644 --- a/PerfectNumbers.Core/DivisorGenerator.cs +++ b/PerfectNumbers.Core/DivisorGenerator.cs @@ -8,31 +8,9 @@ internal static class DivisorGenerator private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool IsValidDivisor( - byte remainder10, - byte remainder8, - byte remainder3, - byte remainder5, - byte remainder7, - byte remainder11, - LastDigit lastDigit) + public static ushort GetDecimalMask(LastDigit lastDigit) { - ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; - if (((decimalMask >> remainder10) & 1) == 0) - { - return false; - } - - if (remainder8 != 1 && remainder8 != 7) - { - return false; - } - - if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0) - { - return false; - } - - return true; + return lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; } + } diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 327d25b..42f64b5 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -36,30 +36,45 @@ internal static class DivisorByDivisorKernels byte firstHitValue; ulong previousExponent = firstExponent; ulong cycleRemainder = 0UL; - if (cycleLength == 0UL) - { - firstHitValue = firstUnity ? (byte)1 : (byte)0; - } - else + bool hasCycle = cycleLength != 0UL; + if (hasCycle) { cycleRemainder = firstExponent % cycleLength; firstHitValue = cycleRemainder == 0UL && firstUnity ? (byte)1 : (byte)0; } + else + { + firstHitValue = firstUnity ? (byte)1 : (byte)0; + } hits[offset] = firstHitValue; byte anyHit = firstHitValue; - for (int i = 1; i < count; i++) + int endIndex = offset + count; + if (hasCycle) { - ulong exponent = exponents[offset + i]; - ulong delta = exponent - previousExponent; - previousExponent = exponent; - - if (cycleLength != 0UL) + for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) { - var remainderValue = new GpuUInt128(cycleRemainder); - remainderValue.AddMod(delta, cycleLength); - cycleRemainder = remainderValue.Low; + ulong exponent = exponents[exponentIndex]; + ulong delta = exponent - previousExponent; + previousExponent = exponent; + + if (delta >= cycleLength) + { + delta -= cycleLength; + if (delta >= cycleLength) + { + delta %= cycleLength; + } + } + + ulong nextRemainder = cycleRemainder + delta; + if (nextRemainder >= cycleLength) + { + nextRemainder -= cycleLength; + } + + cycleRemainder = nextRemainder; byte hit = 0; if (cycleRemainder == 0UL) @@ -67,13 +82,17 @@ internal static class DivisorByDivisorKernels hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; } - hits[offset + i] = hit; + hits[exponentIndex] = hit; anyHit |= hit; } - else + } + else + { + for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) { + ulong exponent = exponents[exponentIndex]; byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; - hits[offset + i] = hit; + hits[exponentIndex] = hit; anyHit |= hit; } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 8d83ee2..5eaf421 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -112,7 +112,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.HitIndexBuffer, resources.Divisors, resources.Exponents, - resources.FilteredDivisors, resources.DivisorData, resources.Offsets, resources.Counts, @@ -177,7 +176,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D hitIndexBuffer, ulong[] divisors, ulong[] exponents, - ulong[] filteredDivisors, GpuDivisorPartialData[] divisorData, int[] offsets, int[] counts, @@ -192,7 +190,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv processedCount = 0UL; lastProcessed = 0UL; - int chunkCapacity = Math.Max(1, batchCapacity); + // EvenPerfectBitScanner always configures non-empty GPU divisor buffers, so the defensive guard stays + // commented out. Tests and benchmarks that rely on zero-capacity buffers must be updated instead of re-enabling + // this branch. + // if (batchCapacity <= 0) + // { + // throw new InvalidOperationException("Divisor buffers must be non-empty."); + // } + + int chunkCountBaseline = batchCapacity; UInt128 twoP128 = (UInt128)prime << 1; UInt128 allowedMax128 = allowedMax; @@ -209,48 +215,63 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; - ulong currentK = 1UL; + ulong remainingCount = maxK; - byte step10 = (byte)(twoP128 % 10UL); - byte step8 = (byte)(twoP128 % 8UL); - byte step5 = (byte)(twoP128 % 5UL); - byte step3 = (byte)(twoP128 % 3UL); - byte step7 = (byte)(twoP128 % 7UL); - byte step11 = (byte)(twoP128 % 11UL); + ulong stepValue = (ulong)twoP128; + byte step10 = (byte)(stepValue % 10UL); + byte step8 = (byte)(stepValue & 7UL); + byte step5 = (byte)(stepValue % 5UL); + byte step3 = (byte)(stepValue % 3UL); + byte step7 = (byte)(stepValue % 7UL); + byte step11 = (byte)(stepValue % 11UL); UInt128 currentDivisor128 = firstDivisor128; - byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL)); - byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL)); - byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL)); - byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL)); - byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL)); - byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL)); + ulong currentDivisorValue = (ulong)currentDivisor128; + byte remainder10 = (byte)(currentDivisorValue % 10UL); + byte remainder8 = (byte)(currentDivisorValue & 7UL); + byte remainder5 = (byte)(currentDivisorValue % 5UL); + 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); - Span filteredStorage = filteredDivisors.AsSpan(); Span divisorStorage = divisors.AsSpan(); Span exponentStorage = exponents.AsSpan(); Span divisorDataStorage = divisorData.AsSpan(); Span offsetStorage = offsets.AsSpan(); Span countStorage = counts.AsSpan(); Span cycleStorage = cycles.AsSpan(); - while (currentK <= maxK && !composite) + Span divisorSpan = divisorStorage; + Span exponentSpan = exponentStorage; + Span divisorDataSpan = divisorDataStorage; + 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; + ArrayView1D countView = countBuffer.View; + ArrayView1D cycleView = cycleBuffer.View; + ArrayView1D exponentViewDevice = exponentBuffer.View; + ArrayView1D hitsView = hitsBuffer.View; + ArrayView1D hitIndexView = hitIndexBuffer.View; + + while (remainingCount > 0UL && !composite) { - int chunkCount = Math.Min(chunkCapacity, batchCapacity); - ulong remainingK = maxK - currentK + 1UL; - if ((ulong)chunkCount > remainingK) + int chunkCount = chunkCountBaseline; + if ((ulong)chunkCount > remainingCount) { - chunkCount = (int)remainingK; + chunkCount = (int)remainingCount; } - if (chunkCount <= 0) - { - processedAll = true; - break; - } - - int filteredCount = 0; UInt128 nextDivisor128 = currentDivisor128; + int admissibleCount = 0; byte localRemainder10 = remainder10; byte localRemainder8 = remainder8; @@ -261,20 +282,32 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv for (int i = 0; i < chunkCount; i++) { - ulong candidate = (ulong)nextDivisor128; - - if (DivisorGenerator.IsValidDivisor(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit)) + bool passesSmallModuli = localRemainder3 != 0 && localRemainder5 != 0 && localRemainder7 != 0 && localRemainder11 != 0; + if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0) { - filteredStorage[filteredCount++] = candidate; + ulong divisorValue = (ulong)nextDivisor128; + MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); + ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData); + if (divisorCycle == prime) + { + int targetIndex = admissibleCount; + divisorSpan[targetIndex] = divisorValue; + exponentSpan[targetIndex] = prime; + divisorDataSpan[targetIndex] = new GpuDivisorPartialData(divisorValue); + offsetSpan[targetIndex] = targetIndex; + countSpan[targetIndex] = 1; + cycleSpan[targetIndex] = divisorCycle; + admissibleCount++; + } } nextDivisor128 += twoP128; - localRemainder10 = AddMod(localRemainder10, step10, 10); - localRemainder8 = AddMod(localRemainder8, step8, 8); - localRemainder5 = AddMod(localRemainder5, step5, 5); - localRemainder3 = AddMod(localRemainder3, step3, 3); - localRemainder7 = AddMod(localRemainder7, step7, 7); - localRemainder11 = AddMod(localRemainder11, step11, 11); + localRemainder10 = AddMod10(localRemainder10, step10); + localRemainder8 = AddMod8(localRemainder8, step8); + localRemainder5 = AddMod5(localRemainder5, step5); + localRemainder3 = AddMod3(localRemainder3, step3); + localRemainder7 = AddMod7(localRemainder7, step7); + localRemainder11 = AddMod11(localRemainder11, step11); } remainder10 = localRemainder10; @@ -289,69 +322,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv UInt128 lastDivisor128 = nextDivisor128 - twoP128; lastProcessed = (ulong)lastDivisor128; currentDivisor128 = nextDivisor128; - - currentK += (ulong)chunkCount; - processedAll = currentK > maxK; - - if (filteredCount == 0) - { - continue; - } - - Span filteredDivisorsSpan = filteredStorage[..filteredCount]; - - Span divisorSpan = divisorStorage; - Span exponentSpan = exponentStorage; - Span divisorDataSpan = divisorDataStorage; - Span offsetSpan = offsetStorage; - Span countSpan = countStorage; - Span cycleSpan = cycleStorage; - int admissibleCount = 0; - int exponentIndex = 0; - for (int i = 0; i < filteredCount; i++) - { - ulong divisorValue = filteredDivisorsSpan[i]; - MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); - ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData); - if (divisorCycle != prime) - { - continue; - } - - divisorSpan[admissibleCount] = divisorValue; - exponentSpan[exponentIndex] = prime; - divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); - offsetSpan[admissibleCount] = exponentIndex; - countSpan[admissibleCount] = 1; - cycleSpan[admissibleCount] = divisorCycle; - - admissibleCount++; - exponentIndex++; - } + remainingCount -= (ulong)chunkCount; if (admissibleCount == 0) { continue; } - ArrayView1D divisorView = divisorDataBuffer.View.SubView(0, admissibleCount); - ArrayView1D offsetView = offsetBuffer.View.SubView(0, admissibleCount); - ArrayView1D countView = countBuffer.View.SubView(0, admissibleCount); - ArrayView1D cycleView = cycleBuffer.View.SubView(0, admissibleCount); - ArrayView1D exponentView = exponentBuffer.View.SubView(0, exponentIndex); - ArrayView1D hitsView = hitsBuffer.View.SubView(0, exponentIndex); - ArrayView1D hitIndexView = hitIndexBuffer.View.SubView(0, 1); - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), admissibleCount); - offsetView.CopyFromCPU(ref MemoryMarshal.GetReference(offsetSpan), admissibleCount); - countView.CopyFromCPU(ref MemoryMarshal.GetReference(countSpan), admissibleCount); - cycleView.CopyFromCPU(ref MemoryMarshal.GetReference(cycleSpan), admissibleCount); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), exponentIndex); + divisorDataView.CopyFromCPU(ref divisorDataRef, admissibleCount); + offsetView.CopyFromCPU(ref offsetRef, admissibleCount); + countView.CopyFromCPU(ref countRef, admissibleCount); + cycleView.CopyFromCPU(ref cycleRef, admissibleCount); + exponentViewDevice.CopyFromCPU(ref exponentRef, admissibleCount); int sentinel = int.MaxValue; hitIndexView.CopyFromCPU(ref sentinel, 1); - kernel(admissibleCount, divisorView, offsetView, countView, exponentView, cycleView, hitsView, hitIndexView); + kernel(admissibleCount, divisorDataView, offsetView, countView, exponentViewDevice, cycleView, hitsView, hitIndexView); int firstHit = sentinel; hitIndexView.CopyToCPU(ref firstHit, 1); @@ -362,6 +349,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; } + processedAll = remainingCount == 0UL; coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; } @@ -379,23 +367,80 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod(byte value, byte delta, byte modulus) + private static byte AddMod3(byte value, byte delta) { - int result = value + delta; - if (result >= modulus) + const int Modulus = 3; + int sum = value + delta; + + if (sum >= Modulus) { - result -= modulus; - if (result >= modulus) - { - result -= modulus; - } + sum -= Modulus; } - return (byte)result; + return (byte)sum; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod5(byte value, byte delta) + { + const int Modulus = 5; + int sum = value + delta; + if (sum >= Modulus) + { + sum -= Modulus; + } + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod7(byte value, byte delta) + { + const int Modulus = 7; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod10(byte value, byte delta) + { + const int Modulus = 10; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod11(byte value, byte delta) + { + const int Modulus = 11; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod8(byte value, byte delta) + { + return (byte)((value + delta) & 7); + } private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) @@ -557,27 +602,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D hitView = _hitBuffer.View; GpuDivisorPartialData partialData = new GpuDivisorPartialData(divisor); - divisorView.SubView(0, 1).CopyFromCPU(ref partialData, 1); + divisorView.CopyFromCPU(ref partialData, 1); int offsetValue = 0; - offsetView.SubView(0, 1).CopyFromCPU(ref offsetValue, 1); + offsetView.CopyFromCPU(ref offsetValue, 1); int countValue = length; - countView.SubView(0, 1).CopyFromCPU(ref countValue, 1); + countView.CopyFromCPU(ref countValue, 1); ulong cycleValue = cycle; - cycleView.SubView(0, 1).CopyFromCPU(ref cycleValue, 1); + cycleView.CopyFromCPU(ref cycleValue, 1); Span hostSpan = _hostBuffer.AsSpan(0, length); primes.CopyTo(hostSpan); + ref ulong hostRef = ref MemoryMarshal.GetReference(hostSpan); - ArrayView1D exponentSlice = exponentView.SubView(0, length); - exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length); - - ArrayView1D hitSliceView = hitView.SubView(0, length); + exponentView.CopyFromCPU(ref hostRef, length); int sentinel = int.MaxValue; - firstHitView.SubView(0, 1).CopyFromCPU(ref sentinel, 1); + firstHitView.CopyFromCPU(ref sentinel, 1); var kernel = _kernel; if (kernel is null) @@ -586,10 +629,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv throw new InvalidOperationException("GPU divisor kernel was not initialized."); } - kernel(1, divisorView, offsetView, countView, exponentSlice, cycleView, hitSliceView, firstHitView); + kernel(1, divisorView, offsetView, countView, exponentView, cycleView, hitView, firstHitView); Span hitSlice = hits.Slice(0, length); - hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), length); + ref byte hitRef = ref MemoryMarshal.GetReference(hitSlice); + hitView.CopyToCPU(ref hitRef, length); Monitor.Exit(_lease.ExecutionLock); } @@ -688,8 +732,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal ulong[] Exponents { get; private set; } = null!; - internal ulong[] FilteredDivisors { get; private set; } = null!; - internal GpuDivisorPartialData[] DivisorData { get; private set; } = null!; internal int[] Offsets { get; private set; } = null!; @@ -764,7 +806,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Divisors = ulongPool.Rent(capacity); Exponents = ulongPool.Rent(capacity); - FilteredDivisors = ulongPool.Rent(capacity); DivisorData = partialPool.Rent(capacity); Offsets = intPool.Rent(capacity); Counts = intPool.Rent(capacity); @@ -803,14 +844,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool partialPool = ThreadStaticPools.GpuDivisorPartialDataPool; ulongPool.Return(Divisors, clearArray: false); ulongPool.Return(Exponents, clearArray: false); - ulongPool.Return(FilteredDivisors, clearArray: false); partialPool.Return(DivisorData, clearArray: false); intPool.Return(Offsets, clearArray: false); intPool.Return(Counts, clearArray: false); ulongPool.Return(Cycles, clearArray: false); Divisors = Array.Empty(); Exponents = Array.Empty(); - FilteredDivisors = Array.Empty(); DivisorData = Array.Empty(); Offsets = Array.Empty(); Counts = Array.Empty(); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 78c60e0..c87a4bc 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -4,6 +4,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using System.Security.Cryptography; using System.Text; +using System.Runtime.InteropServices; using PerfectNumbers.Core.Gpu; using ILGPU; using ILGPU.Runtime; @@ -379,7 +380,8 @@ public class MersenneDivisorCycles for (int i = 0; i < smallLength; i++) { ulong prime = smallPrimes[i]; - if (smallPrimesSquared[i] > remaining) + ulong primeSquared = smallPrimesSquared[i]; + if (primeSquared > remaining) { break; } @@ -438,26 +440,14 @@ public class MersenneDivisorCycles private static void AddFactor(Dictionary counts, ulong factor) { - if (counts.TryGetValue(factor, out int existing)) - { - counts[factor] = existing + 1; - } - else - { - counts[factor] = 1; - } + ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists); + entry = exists ? entry + 1 : 1; } private static void AddFactor(Dictionary counts, ulong factor, int multiplicity) { - if (counts.TryGetValue(factor, out int existing)) - { - counts[factor] = existing + multiplicity; - } - else - { - counts[factor] = multiplicity; - } + ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists); + entry = exists ? entry + multiplicity : multiplicity; } private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) @@ -535,8 +525,8 @@ public class MersenneDivisorCycles private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper) { ulong working = order; - int actual = 0; - while (actual < multiplicity) + int writeIndex = multiplicity; + while (writeIndex > 0) { if (working < prime) { @@ -549,26 +539,20 @@ public class MersenneDivisorCycles break; } - candidateBuffer[actual] = reduced; + writeIndex--; + candidateBuffer[writeIndex] = reduced; working = reduced; - actual++; } + int actual = multiplicity - writeIndex; if (actual == 0) { return; } - Span candidates = candidateBuffer.Slice(0, actual); + Span candidates = candidateBuffer.Slice(writeIndex, actual); Span evaluations = evaluationBuffer.Slice(0, actual); - for (int left = 0, right = actual - 1; left < right; left++, right--) - { - ulong tmp = candidates[left]; - candidates[left] = candidates[right]; - candidates[right] = tmp; - } - stepper.Reset(); evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]); @@ -598,24 +582,28 @@ public class MersenneDivisorCycles // return 2UL; // } + ulong modulus = n; + ulong nMinus1 = modulus - 1UL; + ulong nMinus2 = modulus - 2UL; + while (true) { - ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL; - ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL; + ulong c = (NextRandomUInt64() % nMinus1) + 1UL; + ulong x = (NextRandomUInt64() % nMinus2) + 2UL; ulong y = x; ulong d = 1UL; while (d == 1UL) { - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); + x = AdvancePolynomial(x, c, modulus); + y = AdvancePolynomial(y, c, modulus); + y = AdvancePolynomial(y, c, modulus); ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, n); + d = BinaryGcd(diff, modulus); } - if (d != n) + if (d != modulus) { return d; } diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 09c4435..131df86 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -93,9 +93,15 @@ public static class MersenneNumberDivisorByDivisorTester } else { - while (startIndex > 0 && candidates[startIndex - 1] >= startPrime) + while (startIndex > 0) { - startIndex--; + int previousIndex = startIndex - 1; + if (candidates[previousIndex] < startPrime) + { + break; + } + + startIndex = previousIndex; } } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 9b83a9c..15d6ecc 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -299,22 +299,16 @@ internal static partial class PrimeOrderCalculator Span candidates = candidateBuffer[..actual]; Span evaluations = evaluationBuffer[..actual]; - for (int left = 0, right = actual - 1; left < right; left++, right--) - { - ulong tmp = candidates[left]; - candidates[left] = candidates[right]; - candidates[right] = tmp; - } - stepper.Reset(); - evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]); - for (int j = 1; j < actual; j++) + int last = actual - 1; + evaluations[last] = stepper.InitializeCpuIsUnity(candidates[last]); + for (int j = last - 1; j >= 0; j--) { evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); } - for (int j = actual - 1; j >= 0; j--) + for (int j = 0; j < actual; j++) { if (!evaluations[j]) { @@ -464,16 +458,16 @@ internal static partial class PrimeOrderCalculator } Span candidates = buffer[..actual]; - candidates.Reverse(); stepper.Reset(); - if (stepper.InitializeCpuIsUnity(candidates[0])) + int last = actual - 1; + if (stepper.InitializeCpuIsUnity(candidates[last])) { return true; } - for (int j = 1; j < actual; j++) + for (int j = last - 1; j >= 0; j--) { if (stepper.ComputeNextIsUnity(candidates[j])) { @@ -694,13 +688,16 @@ internal static partial class PrimeOrderCalculator private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) { ulong previous = previousOrder ?? 0UL; - int previousGroup = previousOrder.HasValue ? GetGroup(previousOrder.Value, prime) : 1; bool hasPrevious = previousOrder.HasValue; + ulong threshold1 = prime >> 3; + ulong threshold2 = prime >> 2; + ulong threshold3 = (prime * 3UL) >> 3; + int previousGroup = hasPrevious ? GetGroup(previous, threshold1, threshold2, threshold3) : 1; candidates.Sort((x, y) => { - CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious); - CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious); + CandidateKey keyX = BuildKey(x, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3); + CandidateKey keyY = BuildKey(y, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3); int compare = keyX.Primary.CompareTo(keyY.Primary); if (compare != 0) { @@ -717,9 +714,9 @@ internal static partial class PrimeOrderCalculator }); } - private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious) + private static CandidateKey BuildKey(ulong value, ulong previous, int previousGroup, bool hasPrevious, ulong threshold1, ulong threshold2, ulong threshold3) { - int group = GetGroup(value, prime); + int group = GetGroup(value, threshold1, threshold2, threshold3); if (group == 0) { return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); @@ -769,21 +766,18 @@ internal static partial class PrimeOrderCalculator return groupOffset + (isGe ? 0 : 1); } - private static int GetGroup(ulong value, ulong prime) + private static int GetGroup(ulong value, ulong threshold1, ulong threshold2, ulong threshold3) { - ulong threshold1 = prime >> 3; if (value <= threshold1) { return 1; } - ulong threshold2 = prime >> 2; if (value <= threshold2) { return 2; } - ulong threshold3 = (prime * 3UL) >> 3; if (value <= threshold3) { return 3; @@ -828,9 +822,11 @@ internal static partial class PrimeOrderCalculator } FactorEntry factor = factors[index]; + int factorExponent = factor.Exponent; ulong primeFactor = factor.Value; ulong contribution = 1UL; - for (int exponent = 0; exponent <= factor.Exponent; exponent++) + ulong contributionLimit = factorExponent == 0 ? 0UL : order / primeFactor; + for (int exponent = 0; exponent <= factorExponent; exponent++) { ulong nextDivisor = divisorProduct * contribution; if (nextDivisor > order) @@ -844,12 +840,12 @@ internal static partial class PrimeOrderCalculator return; } - if (exponent == factor.Exponent) + if (exponent == factorExponent) { break; } - if (contribution > order / primeFactor) + if (contribution > contributionLimit) { break; } @@ -973,24 +969,45 @@ internal static partial class PrimeOrderCalculator } Span candidates = buffer[..actual]; - candidates.Reverse(); stepper.Reset(); + int last = actual - 1; - if (powUsed >= powBudget && powBudget > 0) + bool enforceBudget = powBudget > 0; + if (!enforceBudget) + { + powUsed++; + if (stepper.InitializeCpuIsUnity(candidates[last])) + { + return true; + } + + for (int j = last - 1; j >= 0; j--) + { + powUsed++; + if (stepper.ComputeNextIsUnity(candidates[j])) + { + return true; + } + } + + return false; + } + + if (powUsed >= powBudget) { return true; } powUsed++; - if (stepper.InitializeCpuIsUnity(candidates[0])) + if (stepper.InitializeCpuIsUnity(candidates[last])) { return true; } - for (int j = 1; j < actual; j++) + for (int j = last - 1; j >= 0; j--) { - if (powUsed >= powBudget && powBudget > 0) + if (powUsed >= powBudget) { return true; } diff --git a/docs/ByDivisorPathMap.md b/docs/ByDivisorPathMap.md new file mode 100644 index 0000000..85d5a53 --- /dev/null +++ b/docs/ByDivisorPathMap.md @@ -0,0 +1,35 @@ +# --mersenne=bydivisor Path Map + +This temporary map tracks the code paths exercised by the `--mersenne=bydivisor` mode and notes the latest review status for each branch. Entries move from _Queued_ to _Reviewed_ once the surrounding loops have been audited for invariant work. + +## High-level Flow +- `EvenPerfectBitScanner` CLI parsing → _Reviewed_ +- `MersenneNumberDivisorByDivisorTester.Run` → _Reviewed_ +- `EvenPerfectBitScanner.IO.CalculationResultsFile` + - `LoadCandidatesWithinRange` → _Reviewed_ + - `EnumerateCandidates` → _Reviewed_ +- `IMersenneNumberDivisorByDivisorTester` implementations + - `MersenneNumberDivisorByDivisorCpuTester` + - `IsPrime` → _Reviewed_ + - `CheckDivisors` (128-bit path) → _Reviewed_ + - `CheckDivisors64` → _Reviewed_ + - `MersenneNumberDivisorByDivisorGpuTester` + - `IsPrime` → _Reviewed_ + - `CheckDivisors` → _Reviewed_ + - `DivisorByDivisorKernels.CheckKernel` → _Reviewed_ + +- `MersenneDivisorCycles` + - `TryCalculateCycleLengthForExponentCpu` → _Reviewed_ + - `TryFactorIntoCountsInternal` → _Reviewed_ + - `PollardRho64` → _Reviewed_ + - `ReduceOrder` (includes `ProcessReduceOrderPrime`) → _Reviewed_ +- `PrimeOrderCalculator` + - `InitializeStartingOrderCpu` → _Reviewed_ + - `ExponentLoweringCpu` / `ProcessExponentLoweringPrime` → _Reviewed_ + - `ValidateOrderAgainstFactors` / `ValidateOrderForFactor` → _Reviewed_ + - `TryConfirmCandidateCpu` / `CheckCandidateViolation` → _Reviewed_ + - `TryHeuristicFinishCpu` / `BuildCandidates` / `SortCandidates` → _Reviewed_ + +## Next Steps +- Hold for additional guidance on further --mersenne=bydivisor audits; computation branches documented above are fully reviewed. +- Defer result aggregation and reporting reviews until instructed otherwise, keeping the focus on computation-only branches. diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index fd6671e..e2157c9 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -136,7 +136,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic. * [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement. 8. [done] **Refactor by-divisor CPU/GPU testers to reuse heuristics** - * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop while preserving the existing status counters. + * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop and retiring the legacy status counter. * [done] Streamed GPU batches from the same enumerator and Montgomery preparation data, removing the filtered-divisor scratch arrays while continuing to feed the existing kernel and hit bookkeeping. 9. [done] **Route legacy prime checks through heuristics** * [done] Introduce `HeuristicPrimeTester` so heuristic CPU/GPU checks live alongside the legacy `PrimeTester` implementations while callers opt in explicitly. diff --git a/todo_list.txt b/todo_list.txt index 72016e5..9ceb2df 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -41,7 +41,6 @@ Montgomery Primitives Cycle-Aware Modular Arithmetic Status & Progress Reporting -- PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs:296: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly. Cycle Remainder Steppers - PerfectNumbers.Core/CycleRemainderStepper.cs:28: Inline this reset at the call sites so the hot loops reuse struct reinitialization measured fastest in MersenneDivisorCycleLengthGpuBenchmarks, avoiding the extra method call when scanners need to restart stepping. From 968dd89d0d333b2f7916ff4f6cfeb022014aba28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 27 Oct 2025 21:28:41 +0100 Subject: [PATCH 2/2] Remove by-divisor path map document --- docs/ByDivisorPathMap.md | 35 ----------------------------------- 1 file changed, 35 deletions(-) delete mode 100644 docs/ByDivisorPathMap.md diff --git a/docs/ByDivisorPathMap.md b/docs/ByDivisorPathMap.md deleted file mode 100644 index 85d5a53..0000000 --- a/docs/ByDivisorPathMap.md +++ /dev/null @@ -1,35 +0,0 @@ -# --mersenne=bydivisor Path Map - -This temporary map tracks the code paths exercised by the `--mersenne=bydivisor` mode and notes the latest review status for each branch. Entries move from _Queued_ to _Reviewed_ once the surrounding loops have been audited for invariant work. - -## High-level Flow -- `EvenPerfectBitScanner` CLI parsing → _Reviewed_ -- `MersenneNumberDivisorByDivisorTester.Run` → _Reviewed_ -- `EvenPerfectBitScanner.IO.CalculationResultsFile` - - `LoadCandidatesWithinRange` → _Reviewed_ - - `EnumerateCandidates` → _Reviewed_ -- `IMersenneNumberDivisorByDivisorTester` implementations - - `MersenneNumberDivisorByDivisorCpuTester` - - `IsPrime` → _Reviewed_ - - `CheckDivisors` (128-bit path) → _Reviewed_ - - `CheckDivisors64` → _Reviewed_ - - `MersenneNumberDivisorByDivisorGpuTester` - - `IsPrime` → _Reviewed_ - - `CheckDivisors` → _Reviewed_ - - `DivisorByDivisorKernels.CheckKernel` → _Reviewed_ - -- `MersenneDivisorCycles` - - `TryCalculateCycleLengthForExponentCpu` → _Reviewed_ - - `TryFactorIntoCountsInternal` → _Reviewed_ - - `PollardRho64` → _Reviewed_ - - `ReduceOrder` (includes `ProcessReduceOrderPrime`) → _Reviewed_ -- `PrimeOrderCalculator` - - `InitializeStartingOrderCpu` → _Reviewed_ - - `ExponentLoweringCpu` / `ProcessExponentLoweringPrime` → _Reviewed_ - - `ValidateOrderAgainstFactors` / `ValidateOrderForFactor` → _Reviewed_ - - `TryConfirmCandidateCpu` / `CheckCandidateViolation` → _Reviewed_ - - `TryHeuristicFinishCpu` / `BuildCandidates` / `SortCandidates` → _Reviewed_ - -## Next Steps -- Hold for additional guidance on further --mersenne=bydivisor audits; computation branches documented above are fully reviewed. -- Defer result aggregation and reporting reviews until instructed otherwise, keeping the focus on computation-only branches.