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] 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.