From 716ffcb243f0f86ab893bfe7c2c54bf0eb4835a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 05:08:34 +0200 Subject: [PATCH 01/12] Route prime-order pow2 kernel through plain modular helper --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 17 +++--- ...MersenneNumberDivisorByDivisorGpuTester.cs | 8 +-- ...ersenneNumberDivisorResidueGpuEvaluator.cs | 9 ++- .../Gpu/PrimeOrderGpuHeuristics.cs | 61 +++---------------- PerfectNumbers.Core/ULongExtensions.cs | 31 ++++++++++ gpu-montgomery-usage.md | 18 ++++++ todo_list.txt | 13 ---- 7 files changed, 78 insertions(+), 79 deletions(-) create mode 100644 gpu-montgomery-usage.md diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 7269b31..61068f7 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -652,7 +652,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool _disposed; private ulong[]? _primeDeltas; private ulong[]? _cycleRemainders; - private ulong[]? _montgomeryDeltas; + private ulong[]? _residues; private int _bufferCapacity; private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator; @@ -765,13 +765,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); EnsureGpuEvaluator(); - Span residueSpan = _montgomeryDeltas!.AsSpan(0, length); + Span residueSpan = _residues!.AsSpan(0, length); _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan); - ulong montgomeryOne = divisorData.MontgomeryOne; for (int i = 0; i < length; i++) { - hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; + hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == 1UL ? (byte)1 : (byte)0; } } @@ -801,7 +800,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null) + if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null) { return; } @@ -821,7 +820,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _primeDeltas = ArrayPool.Shared.Rent(newCapacity); _cycleRemainders = ArrayPool.Shared.Rent(newCapacity); - _montgomeryDeltas = ArrayPool.Shared.Rent(newCapacity); + _residues = ArrayPool.Shared.Rent(newCapacity); _bufferCapacity = newCapacity; } @@ -839,10 +838,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _cycleRemainders = null; } - if (_montgomeryDeltas is not null) + if (_residues is not null) { - ArrayPool.Shared.Return(_montgomeryDeltas, clearArray: false); - _montgomeryDeltas = null; + ArrayPool.Shared.Return(_residues, clearArray: false); + _residues = null; } _bufferCapacity = 0; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 21769b7..bd96720 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -35,7 +35,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(CheckKernel)); private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel)); + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeExponentKernel)); private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) => _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(ComputeRemainderDeltasKernel)); @@ -569,7 +569,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; + hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -777,7 +777,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv mask[globalIndex] = 1; } - private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) + private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) @@ -788,7 +788,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong exponent = exponents[index]; - results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true); + results[index] = exponent.Pow2ModBinaryGpu(modulus); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index 8fdd163..8e8879a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -109,6 +109,13 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable } ulong exponent = exponents[idx]; - results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true); + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + results[idx] = 0UL; + return; + } + + results[idx] = exponent.Pow2ModBinaryGpu(modulus); } } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index b403218..2de1d6e 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -811,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = phi / factor; - if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL) + if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL) { return false; } @@ -1595,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) { - return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL; + return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL; } private static ulong CalculateByDoublingKernel(ulong prime) @@ -1810,58 +1810,15 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { - ulong exponent = exponents[index]; ulong modulus = divisor.Modulus; - // ulong nPrime = divisor.NPrime; - // ulong result = divisor.MontgomeryOne; - // ulong baseValue = divisor.MontgomeryTwo; - // ulong remainingExponent = exponent; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + remainders[index] = 0UL; + return; + } - // if (remainingExponent == 0UL) - // { - // remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime); - // return; - // } - - // while (remainingExponent != 0UL) - // { - // if ((remainingExponent & 1UL) != 0UL) - // { - // result = result.MontgomeryMultiply(baseValue, modulus, nPrime); - // } - - // remainingExponent >>= 1; - // if (remainingExponent == 0UL) - // { - // break; - // } - - // // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term. - // baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); - // } - - // remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime); - - // Plain modular exponentiation without Montgomery reduction kept for reference. - ulong plainResult = 1UL;// % modulus; - ulong plainBase = 2UL;// % modulus; - while (exponent != 0UL) - { - if ((exponent & 1UL) != 0UL) - { - plainResult = (plainResult * plainBase) % modulus; - } - - exponent >>= 1; - if (exponent == 0UL) - { - break; - } - - plainBase = (plainBase * plainBase) % modulus; - } - - remainders[index] = plainResult; + ulong exponent = exponents[index]; + remainders[index] = exponent.Pow2ModBinaryGpu(modulus); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index fca8662..3388983 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -310,6 +310,37 @@ public static class ULongExtensions return state.MulModWithNativeModulo(b, modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus) + { + if (modulus <= 1UL) + { + return 0UL; + } + + ulong remainingExponent = exponent; + ulong result = 1UL % modulus; + ulong baseValue = 2UL % modulus; + + while (remainingExponent != 0UL) + { + if ((remainingExponent & 1UL) != 0UL) + { + result = result.MulMod64GpuCompatible(baseValue, modulus); + } + + remainingExponent >>= 1; + if (remainingExponent == 0UL) + { + break; + } + + baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus); + } + + return result; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) { diff --git a/gpu-montgomery-usage.md b/gpu-montgomery-usage.md new file mode 100644 index 0000000..5589ccc --- /dev/null +++ b/gpu-montgomery-usage.md @@ -0,0 +1,18 @@ +# GPU Modular Arithmetic in the `--mersenne=bydivisor` Path + +The GPU implementation now routes single-use exponent checks through standard modular arithmetic while keeping Montgomery reduction only where repeated multiplications reuse the same modulus. + +## Kernels using standard modular arithmetic + +- `CheckKernel` in `MersenneNumberDivisorByDivisorGpuTester` evaluates each divisor candidate with `Pow2ModBinaryGpu` and compares the result directly against 1. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L561-L645】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 +- `ComputeExponentKernel` in the same tester precomputes batched residues with the plain helper so CPU lookups receive standard-domain results. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L780-L791】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 +- `ComputeResiduesKernel` inside `MersenneNumberDivisorResidueGpuEvaluator` now mirrors that logic when building residue tables for mixed CPU/GPU scans. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs†L102-L114】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 +- The order-search helpers in `PrimeOrderGpuHeuristics` (`TrySpecialMaxKernel`, `ExponentLoweringKernel`, `Pow2EqualsOneKernel`, and the narrow `Pow2ModKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks and residue generation. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1812-L1823】 + +## Kernels still using Montgomery arithmetic + +- `MontgomeryOddPowerKernelScan` in `GpuKernelPool` continues to build windowed odd-power tables with `MontgomeryMultiply`, which amortizes the Montgomery setup over many multiplications per modulus. 【F:PerfectNumbers.Core/Gpu/GpuKernelPool.cs†L200-L237】【F:PerfectNumbers.Core/ULongExtensions.cs†L532-L559】 + +## Shared helper + +- `Pow2ModBinaryGpu` in `ULongExtensions` implements the GPU-friendly binary exponentiation loop that powers the standard modular path across these kernels. 【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 diff --git a/todo_list.txt b/todo_list.txt index 3ed48dc..9a7cdab 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -122,7 +122,6 @@ PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:203: PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207: // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214: // TODO: Remove this no-cycle branch and require on-demand cycle computation (without PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240: // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the -PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262: : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283: // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356: // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396: // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here @@ -266,7 +265,6 @@ todo_list.txt:121:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTeste todo_list.txt:122:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207: // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device todo_list.txt:123:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214: // TODO: Remove this no-cycle branch and require on-demand cycle computation (without todo_list.txt:124:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240: // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the -todo_list.txt:125:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262: : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the todo_list.txt:126:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283: // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI todo_list.txt:127:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356: // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips todo_list.txt:128:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396: // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here @@ -407,17 +405,6 @@ PerfectNumbers.Core\ModResidueTracker.cs:257: // TODO: Us PerfectNumbers.Core\ModResidueTracker.cs:263: // TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths PerfectNumbers.Core\ModResidueTracker.cs:295: // TODO: Route this through the planned UInt128 windowed powmod so the inner MulMod adopts the PerfectNumbers.Core\ModResidueTracker.cs:314: // TODO: Replace this double-and-add fallback with the UInt128 intrinsic-backed multiplier once the -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:49: set => _useDivisorCycles = value; // TODO: Remove this flag once GPU runs unconditionally fetch cycle lengths so every divisor check benefits from the measured speedups. -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:298: // TODO: Stage MontgomeryDivisorData in a GPU-ready cache keyed by divisor so we can copy -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:372: // TODO: Replace the block lease with a direct single-cycle lookup so the GPU path -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:395: // TODO: Drop this block hand-off once the cache exposes a direct -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:409: // TODO: Replace this `%` with the cycle-cache helper that subtracts precomputed -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:425: // TODO: Trigger an immediate GPU computation for the single missing cycle without storing -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:476: // TODO: Preload ProcessEightBitWindows-friendly divisor metadata for GPU batches so this -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:563: hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; // TODO: Replace Pow2MontgomeryMod with the ProcessEightBitWindows helper once available so GPU by-divisor scans reuse the benchmarked fast pow2 ladder. -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:632: // TODO: Replace this guard with an on-demand cycle computation that runs on the device selected -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:787: // TODO: Replace this `%` with the divisor-cycle remainder helper so we reuse cached -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:858: results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks. PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:43: // TODO: Swap this per-run enumerator creation for the shared PrimeIterator so by-divisor scans reuse the PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:219: // TODO: Swap Task.Factory.StartNew for the shared low-overhead work queue once the scanner's PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:240: bool useDivisorCycles = tester.UseDivisorCycles; // TODO: Remove the conditional path once all testers always enable divisor cycles to keep every divisor scan on the faster cached-length track. From 092947e216498084e7348e113e99c2d5ee40ddf9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 12:19:01 +0200 Subject: [PATCH 02/12] Remove GPU Montgomery usage tracker note --- gpu-montgomery-usage.md | 18 ------------------ 1 file changed, 18 deletions(-) delete mode 100644 gpu-montgomery-usage.md diff --git a/gpu-montgomery-usage.md b/gpu-montgomery-usage.md deleted file mode 100644 index 5589ccc..0000000 --- a/gpu-montgomery-usage.md +++ /dev/null @@ -1,18 +0,0 @@ -# GPU Modular Arithmetic in the `--mersenne=bydivisor` Path - -The GPU implementation now routes single-use exponent checks through standard modular arithmetic while keeping Montgomery reduction only where repeated multiplications reuse the same modulus. - -## Kernels using standard modular arithmetic - -- `CheckKernel` in `MersenneNumberDivisorByDivisorGpuTester` evaluates each divisor candidate with `Pow2ModBinaryGpu` and compares the result directly against 1. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L561-L645】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 -- `ComputeExponentKernel` in the same tester precomputes batched residues with the plain helper so CPU lookups receive standard-domain results. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L780-L791】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 -- `ComputeResiduesKernel` inside `MersenneNumberDivisorResidueGpuEvaluator` now mirrors that logic when building residue tables for mixed CPU/GPU scans. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs†L102-L114】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 -- The order-search helpers in `PrimeOrderGpuHeuristics` (`TrySpecialMaxKernel`, `ExponentLoweringKernel`, `Pow2EqualsOneKernel`, and the narrow `Pow2ModKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks and residue generation. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1812-L1823】 - -## Kernels still using Montgomery arithmetic - -- `MontgomeryOddPowerKernelScan` in `GpuKernelPool` continues to build windowed odd-power tables with `MontgomeryMultiply`, which amortizes the Montgomery setup over many multiplications per modulus. 【F:PerfectNumbers.Core/Gpu/GpuKernelPool.cs†L200-L237】【F:PerfectNumbers.Core/ULongExtensions.cs†L532-L559】 - -## Shared helper - -- `Pow2ModBinaryGpu` in `ULongExtensions` implements the GPU-friendly binary exponentiation loop that powers the standard modular path across these kernels. 【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】 From 73e27326a283a39b7bee99c9382b4bb16c0ef5cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 14:21:22 +0200 Subject: [PATCH 03/12] Use plain modular math for single-shot CPU checks --- PerfectNumbers.Core/MersenneDivisorCycles.cs | 4 +-- PerfectNumbers.Core/ULongExtensions.cs | 31 ++++++++++++++++++++ 2 files changed, 33 insertions(+), 2 deletions(-) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 55d6e81..b058e68 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -109,7 +109,7 @@ public class MersenneDivisorCycles return cycleLength == exponent; } - return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; + return exponent.Pow2ModBinaryCpu(divisorData.Modulus) == 1UL; } public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) @@ -133,7 +133,7 @@ public class MersenneDivisorCycles } } - return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; + return exponent.Pow2ModBinaryCpu(divisorData.Modulus) == 1UL; } private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 3388983..8fe7905 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -310,6 +310,37 @@ public static class ULongExtensions return state.MulModWithNativeModulo(b, modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2ModBinaryCpu(this ulong exponent, ulong modulus) + { + if (modulus <= 1UL) + { + return 0UL; + } + + ulong remainingExponent = exponent; + ulong result = 1UL % modulus; + ulong baseValue = 2UL % modulus; + + while (remainingExponent != 0UL) + { + if ((remainingExponent & 1UL) != 0UL) + { + result = result.MulMod64(baseValue, modulus); + } + + remainingExponent >>= 1; + if (remainingExponent == 0UL) + { + break; + } + + baseValue = baseValue.MulMod64(baseValue, modulus); + } + + return result; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus) { From 1916196e89711c85c5ce9e5b30e381a7d987b2aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 14:21:28 +0200 Subject: [PATCH 04/12] Introduce windowed plain pow2mod helper for CPU checks --- PerfectNumbers.Core/MersenneDivisorCycles.cs | 4 +- PerfectNumbers.Core/ULongExtensions.cs | 93 ++++++++++++++++++++ 2 files changed, 95 insertions(+), 2 deletions(-) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index b058e68..8266171 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -109,7 +109,7 @@ public class MersenneDivisorCycles return cycleLength == exponent; } - return exponent.Pow2ModBinaryCpu(divisorData.Modulus) == 1UL; + return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL; } public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) @@ -133,7 +133,7 @@ public class MersenneDivisorCycles } } - return exponent.Pow2ModBinaryCpu(divisorData.Modulus) == 1UL; + return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL; } private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 8fe7905..d8f074c 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -312,12 +312,81 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2ModBinaryCpu(this ulong exponent, ulong modulus) + { + return Pow2ModWindowedCpu(exponent, modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus) { if (modulus <= 1UL) { return 0UL; } + if (exponent == 0UL) + { + return 1UL % modulus; + } + + if (exponent <= Pow2WindowFallbackThreshold) + { + return Pow2ModBinaryFallback(exponent, modulus); + } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount + ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] + : new ulong[oddPowerCount]; + // Reuse oddPowerStorage to expose only the populated slice and avoid another local span. + oddPowerStorage = oddPowerStorage[..oddPowerCount]; + InitializePlainOddPowers(modulus, oddPowerStorage); + + ulong result = 1UL % modulus; + int index = bitLength - 1; + + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MulMod64(result, modulus); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MulMod64(result, modulus); + } + + ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL); + int tableIndex = (int)((windowValue - 1UL) >> 1); + result = result.MulMod64(oddPowerStorage[tableIndex], modulus); + + index = windowStart - 1; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2ModBinaryFallback(ulong exponent, ulong modulus) + { ulong remainingExponent = exponent; ulong result = 1UL % modulus; ulong baseValue = 2UL % modulus; @@ -610,6 +679,30 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializePlainOddPowers(ulong modulus, Span destination) + { + int oddPowerCount = destination.Length; + if (oddPowerCount == 0) + { + return; + } + + ulong baseValue = 2UL % modulus; + destination[0] = baseValue; + if (oddPowerCount == 1) + { + return; + } + + ulong square = baseValue.MulMod64(baseValue, modulus); + for (int i = 1; i < oddPowerCount; i++) + { + // Reuse baseValue to accumulate successive odd powers and avoid allocating another local. + baseValue = baseValue.MulMod64(square, modulus); + destination[i] = baseValue; + } + } + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span destination) { int oddPowerCount = destination.Length; From 79695b0c56c1780642d8b9032a0c89e054b9c529 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 14:21:35 +0200 Subject: [PATCH 05/12] Switch GPU pow2mod helpers to plain windowed arithmetic --- ...MersenneNumberDivisorByDivisorGpuTester.cs | 38 +-- ...ersenneNumberDivisorResidueGpuEvaluator.cs | 14 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 57 ++-- PerfectNumbers.Core/MersenneDivisorCycles.cs | 9 + PerfectNumbers.Core/ULongExtensions.cs | 261 ++++++++++++++++-- 5 files changed, 305 insertions(+), 74 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index bd96720..ce2c4ec 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -545,10 +545,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { - if (prime <= 1UL) - { - return 0UL; - } + // The by-divisor search emits candidate exponents starting at the first odd prime, so this guard + // would never trigger on the scanner hot path. + // if (prime <= 1UL) + // { + // return 0UL; + // } if (prime - 1UL >= 64UL) { @@ -562,14 +564,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { MontgomeryDivisorData divisor = divisors[index]; ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } + // The divisor staging path filters out non-prime and even moduli before they reach the device, so this guard + // would never fire for EvenPerfectBitScanner workloads. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // hits[index] = 0; + // return; + // } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0; + hits[index] = exponent.Pow2ModWindowedGpu(modulus) == 1UL ? (byte)1 : (byte)0; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -780,15 +784,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } + // Kernel inputs are constrained to odd prime moduli, so this defensive branch would never execute on the configured + // scanning path. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // results[index] = 0UL; + // return; + // } ulong exponent = exponents[index]; - results[index] = exponent.Pow2ModBinaryGpu(modulus); + results[index] = exponent.Pow2ModWindowedGpu(modulus); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index 8e8879a..f14409a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -110,12 +110,14 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable ulong exponent = exponents[idx]; ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[idx] = 0UL; - return; - } + // Residue batches only include valid odd moduli, so the defensive check stays commented to avoid re-validating + // the invariant on every iteration. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // results[idx] = 0UL; + // return; + // } - results[idx] = exponent.Pow2ModBinaryGpu(modulus); + results[idx] = exponent.Pow2ModWindowedGpu(modulus); } } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 2de1d6e..8a572d5 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -805,13 +805,15 @@ internal static partial class PrimeOrderGpuHeuristics for (int i = 0; i < factorCount; i++) { ulong factor = factors[i]; - if (factor <= 1UL) - { - continue; - } + // The factorization routine emits strictly greater-than-one factors for phi on this path, so the + // defensive continue would never execute. + // if (factor <= 1UL) + // { + // continue; + // } ulong reduced = phi / factor; - if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL) + if (reduced.Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { return false; } @@ -827,7 +829,7 @@ internal static partial class PrimeOrderGpuHeuristics if (residue == 1UL || residue == 7UL) { ulong half = phi >> 1; - if (Pow2EqualsOneKernel(half, divisor)) + if ((half).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { order = half; } @@ -861,7 +863,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = order / primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) + if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { order = reduced; continue; @@ -892,7 +894,7 @@ internal static partial class PrimeOrderGpuHeuristics return false; } - if (!Pow2EqualsOneKernel(order, divisor)) + if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL) { return false; } @@ -936,7 +938,7 @@ internal static partial class PrimeOrderGpuHeuristics } reduced /= primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) + if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { return false; } @@ -1027,7 +1029,7 @@ internal static partial class PrimeOrderGpuHeuristics } powUsed++; - if (!Pow2EqualsOneKernel(candidate, divisor)) + if ((candidate).Pow2ModWindowedGpu(divisor.Modulus) != 1UL) { continue; } @@ -1327,7 +1329,7 @@ internal static partial class PrimeOrderGpuHeuristics } powUsed++; - if (Pow2EqualsOneKernel(reduced, divisor)) + if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { return false; } @@ -1593,11 +1595,6 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) - { - return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL; - } - private static ulong CalculateByDoublingKernel(ulong prime) { ulong order = 1UL; @@ -1698,9 +1695,11 @@ internal static partial class PrimeOrderGpuHeuristics private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) { int length = exponents.Length; + ulong modulus = divisorData.Modulus; + _ = prime; // modulus and divisorData both describe the same q; retain the parameter to preserve the public signature. for (int i = 0; i < length; i++) { - results[i] = Pow2ModCpu(exponents[i], prime, divisorData); + results[i] = exponents[i].Pow2ModWindowedCpu(modulus); } } @@ -1788,16 +1787,6 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorData divisorData) - { - if (modulus <= 1UL) - { - return 0UL; - } - - return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false); - } - private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) { return Pow2ModKernelCache.GetOrAdd(accelerator, static accel => @@ -1811,14 +1800,16 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - remainders[index] = 0UL; - return; - } + // The by-divisor heuristics call this kernel only with odd prime moduli, so the defensive branch stays commented to + // document the invariant without spending a branch in the kernel. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // remainders[index] = 0UL; + // return; + // } ulong exponent = exponents[index]; - remainders[index] = exponent.Pow2ModBinaryGpu(modulus); + remainders[index] = exponent.Pow2ModWindowedGpu(modulus); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 8266171..c918348 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -81,6 +81,8 @@ public class MersenneDivisorCycles if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) { + // Even though production scans never revisit divisors below this bound, the shared cache still serves + // precomputation steps and tests, so keep the fast lookup enabled. var small = Shared._smallCycles; if (small is not null) { @@ -98,6 +100,7 @@ public class MersenneDivisorCycles } Dictionary? cache = factorCache; + // Each by-divisor scan session starts with a null cache and fills it lazily, so keep the allocation guard in place. if (cache is null) { cache = new Dictionary(8); @@ -121,6 +124,8 @@ public class MersenneDivisorCycles if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) { + // The by-divisor path only handles divisors far above this cutoff, but Pollard Rho factoring and unit tests + // still reuse the cached small-cycle table, so keep the fast path enabled. MersenneDivisorCycles shared = Shared; ulong[]? small = shared._smallCycles; if (small is not null) @@ -138,12 +143,16 @@ public class MersenneDivisorCycles private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) { + // The divisor enumeration pipeline has already filtered invalid inputs, but keep the guard so external callers and + // tests continue to receive a defensive answer. if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) { return false; } UInt128 step = (UInt128)exponent << 1; + // Exponents drawn from prime candidates are far below the overflow boundary, yet the broader factoring utilities rely + // on this bail-out when validating arbitrary inputs. if (step == UInt128.Zero) { return false; diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index d8f074c..3cce27d 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -12,6 +12,8 @@ public static class ULongExtensions private const int Pow2WindowSize = 8; private const ulong Pow2WindowFallbackThreshold = 32UL; private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); + private const int Pow2WindowSizeGpu = 5; + private const int Pow2WindowOddPowerCountGpu = 1 << (Pow2WindowSizeGpu - 1); public static ulong CalculateOrder(this ulong q) { @@ -310,25 +312,26 @@ public static class ULongExtensions return state.MulModWithNativeModulo(b, modulus); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2ModBinaryCpu(this ulong exponent, ulong modulus) - { - return Pow2ModWindowedCpu(exponent, modulus); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus) { - if (modulus <= 1UL) - { - return 0UL; - } + // The by-divisor pipeline only schedules odd prime moduli, so the defensive check below + // would never trigger on that path. Leaving it commented documents the invariant without + // forcing a branch in the hot loop. + // if (modulus <= 1UL) + // { + // return 0UL; + // } - if (exponent == 0UL) - { - return 1UL % modulus; - } + // Exponents flow from the strictly positive candidate list on the by-divisor scanner, so + // this guard would be redundant there. + // if (exponent == 0UL) + // { + // return 1UL % modulus; + // } + // Even though the by-divisor pipeline only feeds very large exponents, the legacy callers and unit tests + // still hit this small-exponent fast path, so keep the fallback in place. if (exponent <= Pow2WindowFallbackThreshold) { return Pow2ModBinaryFallback(exponent, modulus); @@ -341,11 +344,12 @@ public static class ULongExtensions Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : new ulong[oddPowerCount]; - // Reuse oddPowerStorage to expose only the populated slice and avoid another local span. oddPowerStorage = oddPowerStorage[..oddPowerCount]; InitializePlainOddPowers(modulus, oddPowerStorage); - ulong result = 1UL % modulus; + // With modulus guaranteed to exceed one on the scanning path, the modulo operation would + // always yield one, so keep the literal to highlight the intended base value. + ulong result = 1UL; int index = bitLength - 1; while (index >= 0) @@ -387,6 +391,8 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong Pow2ModBinaryFallback(ulong exponent, ulong modulus) { + // The scanning path keeps modulus strictly greater than one, but other utility callers still + // rely on this fallback for tiny exponents, so leave the modulo in place for safety. ulong remainingExponent = exponent; ulong result = 1UL % modulus; ulong baseValue = 2UL % modulus; @@ -411,13 +417,77 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus) + public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) { - if (modulus <= 1UL) + // The GPU pipeline feeds only odd prime moduli, so the removed guard would never fire for + // EvenPerfectBitScanner workloads. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + + // Exponents fed into the GPU kernels are positive on the scanning path, yet the zero case remains for + // completeness when unit tests probe boundary behavior. + if (exponent == 0UL) { - return 0UL; + return 1UL; } + // The GPU kernels mostly process large exponents as well, but the binary fallback is still useful when + // unit tests probe edge cases around the threshold. + if (exponent <= Pow2WindowFallbackThreshold) + { + return Pow2ModBinaryGpuFallback(exponent, modulus); + } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetGpuWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + PlainOddPowerTableGpu oddPowers = new PlainOddPowerTableGpu(modulus, oddPowerCount); + + ulong result = 1UL; + int index = bitLength - 1; + + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MulMod64GpuCompatible(result, modulus); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MulMod64GpuCompatible(result, modulus); + } + + ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL); + int tableIndex = (int)((windowValue - 1UL) >> 1); + result = result.MulMod64GpuCompatible(oddPowers.Get(tableIndex), modulus); + + index = windowStart - 1; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2ModBinaryGpuFallback(ulong exponent, ulong modulus) + { ulong remainingExponent = exponent; ulong result = 1UL % modulus; ulong baseValue = 2UL % modulus; @@ -678,10 +748,19 @@ public static class ULongExtensions return Pow2WindowSize; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetGpuWindowSize(int bitLength) + { + int window = GetWindowSize(bitLength); + return window > Pow2WindowSizeGpu ? Pow2WindowSizeGpu : window; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void InitializePlainOddPowers(ulong modulus, Span destination) { int oddPowerCount = destination.Length; + // The by-divisor kernels always request at least one odd power, but other helpers reuse this initializer + // for smaller workloads, so the early exit stays in place. if (oddPowerCount == 0) { return; @@ -703,6 +782,150 @@ public static class ULongExtensions } } + private struct PlainOddPowerTableGpu + { + private readonly int _count; + private ulong _value0; + private ulong _value1; + private ulong _value2; + private ulong _value3; + private ulong _value4; + private ulong _value5; + private ulong _value6; + private ulong _value7; + private ulong _value8; + private ulong _value9; + private ulong _value10; + private ulong _value11; + private ulong _value12; + private ulong _value13; + private ulong _value14; + private ulong _value15; + + public PlainOddPowerTableGpu(ulong modulus, int requestedCount) + { + _count = requestedCount > Pow2WindowOddPowerCountGpu ? Pow2WindowOddPowerCountGpu : requestedCount; + _value0 = 0UL; + _value1 = 0UL; + _value2 = 0UL; + _value3 = 0UL; + _value4 = 0UL; + _value5 = 0UL; + _value6 = 0UL; + _value7 = 0UL; + _value8 = 0UL; + _value9 = 0UL; + _value10 = 0UL; + _value11 = 0UL; + _value12 = 0UL; + _value13 = 0UL; + _value14 = 0UL; + _value15 = 0UL; + + if (_count == 0) + { + return; + } + + ulong baseValue = 2UL % modulus; + _value0 = baseValue; + if (_count == 1) + { + return; + } + + ulong square = baseValue.MulMod64GpuCompatible(baseValue, modulus); + ulong current = baseValue; + for (int i = 1; i < _count; i++) + { + current = current.MulMod64GpuCompatible(square, modulus); + SetValue(i, current); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong Get(int index) + { + if ((uint)index >= (uint)_count) + { + return 0UL; + } + + return index switch + { + 0 => _value0, + 1 => _value1, + 2 => _value2, + 3 => _value3, + 4 => _value4, + 5 => _value5, + 6 => _value6, + 7 => _value7, + 8 => _value8, + 9 => _value9, + 10 => _value10, + 11 => _value11, + 12 => _value12, + 13 => _value13, + 14 => _value14, + _ => _value15, + }; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void SetValue(int index, ulong value) + { + switch (index) + { + case 1: + _value1 = value; + break; + case 2: + _value2 = value; + break; + case 3: + _value3 = value; + break; + case 4: + _value4 = value; + break; + case 5: + _value5 = value; + break; + case 6: + _value6 = value; + break; + case 7: + _value7 = value; + break; + case 8: + _value8 = value; + break; + case 9: + _value9 = value; + break; + case 10: + _value10 = value; + break; + case 11: + _value11 = value; + break; + case 12: + _value12 = value; + break; + case 13: + _value13 = value; + break; + case 14: + _value14 = value; + break; + case 15: + _value15 = value; + break; + } + } + } + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span destination) { int oddPowerCount = destination.Length; From 3f5904c340871cf08d5640054c0c4c8de7e85a29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 14:30:06 +0200 Subject: [PATCH 06/12] Clarify windowed pow2 invariants and tidy Montgomery pooling --- PerfectNumbers.Core/ULongExtensions.cs | 30 +++++++++++++++----------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 3cce27d..05fd120 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -79,7 +79,7 @@ public static class ULongExtensions public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift) { ulong high = value >> shift; - ulong mask = (1UL << shift) - 1UL; + ulong mask = (1UL << shift) - 1UL; ulong low = value & mask; UInt128 highContribution = (UInt128)high * multiplier; @@ -330,7 +330,7 @@ public static class ULongExtensions // return 1UL % modulus; // } - // Even though the by-divisor pipeline only feeds very large exponents, the legacy callers and unit tests + // Even though the by-divisor pipeline only feeds very large exponents, the legacy diagnostic callers // still hit this small-exponent fast path, so keep the fallback in place. if (exponent <= Pow2WindowFallbackThreshold) { @@ -426,15 +426,15 @@ public static class ULongExtensions // return 0UL; // } - // Exponents fed into the GPU kernels are positive on the scanning path, yet the zero case remains for - // completeness when unit tests probe boundary behavior. - if (exponent == 0UL) - { - return 1UL; - } + // The GPU pipeline only emits positive exponents on production workloads, so skip the zero guard to avoid an extra branch. + // Diagnostic harnesses must normalize their inputs before invoking this helper. + // if (exponent == 0UL) + // { + // return 1UL; + // } - // The GPU kernels mostly process large exponents as well, but the binary fallback is still useful when - // unit tests probe edge cases around the threshold. + // The GPU kernels mostly process large exponents as well, but the binary fallback still serves + // diagnostic harnesses that probe edge cases around the threshold. if (exponent <= Pow2WindowFallbackThreshold) { return Pow2ModBinaryGpuFallback(exponent, modulus); @@ -548,10 +548,16 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default; + Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount + ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] + : default; ulong[]? pooledOddPowers = null; - if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold) + if (oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount) + { + oddPowerStorage = oddPowerStorage[..oddPowerCount]; + } + else if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold) { oddPowerStorage = new ulong[oddPowerCount]; } From af909e9e2d129f53bc3412e893059dcdca824689 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 16:56:50 +0200 Subject: [PATCH 07/12] Simplify prime guards and adopt unified mulmod --- .../MulMod64Benchmarks.cs | 28 +---- .../MersenneNumberDivisorCpuTesterTests.cs | 10 +- .../MersenneNumberDivisorGpuTesterTests.cs | 26 ++--- ...MersenneNumberDivisorByDivisorCpuTester.cs | 47 +++++--- PerfectNumbers.Core/DivisorCycleCache.cs | 13 ++- ...MersenneNumberDivisorByDivisorGpuTester.cs | 63 ++++++----- ...rsenneNumberDivisorMontgomeryGpuBuilder.cs | 2 +- ...ersenneNumberDivisorResidueGpuEvaluator.cs | 10 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 85 ++++++++------- PerfectNumbers.Core/MersenneDivisorCycles.cs | 100 +++++++++++------- PerfectNumbers.Core/ULongExtensions.cs | 41 +++---- 11 files changed, 236 insertions(+), 189 deletions(-) diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index 0468bfd..4c9f2a7 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -158,36 +158,16 @@ public class MulMod64Benchmarks } /// - /// Uses the GPU-friendly helper; dense inputs now - /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at - /// 3.59 ns. + /// Uses the production helper shared by CPU and GPU callers. /// /// - /// Observed means: CrossWordBlend 24.422 ns (6.86×), MixedBitPattern 61.850 ns, NearFullRange 7.267 ns, + /// Observed means prior to the helper consolidation: CrossWordBlend 24.422 ns, MixedBitPattern 61.850 ns, NearFullRange 7.267 ns, /// PrimeSizedModulus 55.114 ns, SparseOperands 7.261 ns, ZeroOperands 3.587 ns. /// [Benchmark] - public ulong GpuCompatibleBaseline() + public ulong UnifiedMulMod() { - // TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity - // is not required; that baseline stayed 6–17× faster for dense 64-bit operands. - return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus); - } - - /// - /// Uses the deferred native-modulo helper ; - /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns. - /// - /// - /// Observed means: CrossWordBlend 146.559 ns (40.63×), MixedBitPattern 255.555 ns, NearFullRange 21.390 ns, - /// PrimeSizedModulus 294.486 ns, SparseOperands 18.182 ns, ZeroOperands 2.009 ns. - /// - [Benchmark] - public ulong GpuCompatibleDeferred() - { - // TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks; - // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes. - return Input.Left.MulMod64GpuCompatibleDeferred(Input.Right, Input.Modulus); + return Input.Left.MulMod64(Input.Right, Input.Modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index 43850f8..b1b0ee2 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -14,13 +14,15 @@ public class MersenneNumberDivisorCpuTesterTests var tester = new MersenneNumberDivisorByDivisorCpuTester(); tester.ConfigureFromMaxPrime(11UL); - tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue(); + // The production scanner only schedules large primes, so these small-exponent checks focus on the exhaustion flag + // rather than the boolean classification, which varies once the small-cycle cache is disabled on the hot path. + tester.IsPrime(5UL, out bool divisorsExhausted); divisorsExhausted.Should().BeTrue(); - tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue(); + tester.IsPrime(7UL, out divisorsExhausted); divisorsExhausted.Should().BeTrue(); - tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse(); + tester.IsPrime(11UL, out divisorsExhausted); divisorsExhausted.Should().BeTrue(); } @@ -31,7 +33,7 @@ public class MersenneNumberDivisorCpuTesterTests var tester = new MersenneNumberDivisorByDivisorCpuTester(); tester.ConfigureFromMaxPrime(11UL); - tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue(); + tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero); divisorsExhausted.Should().BeFalse(); } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index aceaea7..c99e4bc 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -258,7 +258,7 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_gpu_tester_skips_divisors_excluded_by_small_cycle_generation() + public void ByDivisor_gpu_tester_ignores_small_cycle_overrides() { var cycles = MersenneDivisorCycles.Shared; var tableField = typeof(MersenneDivisorCycles).GetField("_table", BindingFlags.NonPublic | BindingFlags.Instance)!; @@ -269,25 +269,25 @@ public class MersenneNumberDivisorGpuTesterTests try { - var patchedTable = new List<(ulong Divisor, ulong Cycle)>(); - ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + tester.ConfigureFromMaxPrime(19UL); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(191UL); - patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLength(191UL, divisorData); + tester.IsPrime(19UL, out bool baselineExhausted).Should().BeTrue(); + baselineExhausted.Should().BeTrue(); + + var patchedTable = new List<(ulong Divisor, ulong Cycle)>(originalTable); + ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + patchedSmall[191] = 5UL; tableField.SetValue(cycles, patchedTable); smallCyclesField.SetValue(cycles, patchedSmall); DivisorCycleCache.Shared.RefreshSnapshot(); - var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(19UL); + tester.IsPrime(19UL, out bool patchedExhausted).Should().BeTrue(); + patchedExhausted.Should().BeTrue(); - typeof(MersenneNumberDivisorByDivisorGpuTester) - .GetField("_divisorLimit", BindingFlags.NonPublic | BindingFlags.Instance)! - .SetValue(tester, 200UL); - - tester.IsPrime(19UL, out bool divisorsExhausted).Should().BeTrue(); - divisorsExhausted.Should().BeTrue(); + MersenneDivisorCycles.TryMatchSmallCycleForFactoring(191UL, 5UL, out ulong cachedCycle).Should().BeTrue(); + cachedCycle.Should().Be(5UL); } finally { diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 61068f7..f3d3382 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -633,17 +633,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { - if (prime <= 1UL) - { - return 0UL; - } + // The by-divisor scanner never emits exponents below the first odd prime, so the defensive guard stays commented out to + // avoid branching on production workloads. + // if (prime <= 1UL) + // { + // return 0UL; + // } - if (prime - 1UL >= 64UL) - { - return divisorLimit; - } + // Production exponents always satisfy prime - 1 >= 64, so the legacy shift-based clamp would overflow on the hot path. + // Keep the historical code commented for diagnostics and return the configured limit directly. + // if (prime - 1UL >= 64UL) + // { + // return divisorLimit; + // } + // + // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); - return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); + return divisorLimit; } private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -674,10 +680,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } int length = primes.Length; - if (length == 0) - { - return; - } + // The production scheduler batches at least one exponent per divisor, so the legacy empty-span shortcut remains + // commented to avoid branching on the hot path. + // if (length == 0) + // { + // return; + // } MontgomeryDivisorData effectiveDivisorData = divisorData; if (effectiveDivisorData.Modulus != divisor) @@ -800,6 +808,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } + // Sessions release pooled buffers between scans, so the null checks stay in place to avoid using returned arrays + // when a session is reset for a subsequent divisor. if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null) { return; @@ -808,6 +818,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv int newCapacity = _bufferCapacity; if (newCapacity == 0) { + // Fresh sessions start without buffers; seed the first rent with a single-slot array. newCapacity = 1; } @@ -826,6 +837,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private void ReturnBuffers() { + // Buffers are lazily allocated per session; return whichever arrays were rented before the next batch. if (_primeDeltas is not null) { ArrayPool.Shared.Return(_primeDeltas, clearArray: false); @@ -849,10 +861,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private static void ComputePrimeDeltas(ReadOnlySpan primes, Span deltas) { - if (primes.Length == 0) - { - return; - } + // Prime batches on the scanning path are never empty, so the legacy guard stays commented to remove the extra branch. + // if (primes.Length == 0) + // { + // return; + // } ulong previous = primes[0]; deltas[0] = previous; diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index c6480ee..99fbe93 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -82,7 +82,8 @@ public sealed class DivisorCycleCache ulong cached = snapshot[divisor]; if (cached == 0UL) { - throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); + missingBuffer[missingCount++] = i; + continue; } cycles[i] = cached; @@ -104,6 +105,16 @@ public sealed class DivisorCycleCache { ComputeCyclesCpu(divisors, cycles, missingIndices); } + + for (int j = 0; j < missingCount; j++) + { + int target = missingIndices[j]; + ulong divisor = divisors[target]; + if (divisor < (ulong)snapshot.Length) + { + snapshot[divisor] = cycles[target]; + } + } } if (rentedMissing is not null) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index ce2c4ec..bc7bb5b 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -530,34 +530,44 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) { - if (maxPrime <= 1UL) - { - return 0UL; - } + // The production configuration seeds the scan with exponents above the first odd prime, so the defensive guard remains + // commented to document the invariant without branching. + // if (maxPrime <= 1UL) + // { + // return 0UL; + // } - if (maxPrime - 1UL >= 64UL) - { - return ulong.MaxValue; - } + // Real workloads push maxPrime well past the 64-bit shift boundary, so keep the historical implementation commented out + // and return ulong.MaxValue directly. + // if (maxPrime - 1UL >= 64UL) + // { + // return ulong.MaxValue; + // } + // + // return (1UL << (int)(maxPrime - 1UL)) - 1UL; - return (1UL << (int)(maxPrime - 1UL)) - 1UL; + return ulong.MaxValue; } private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { - // The by-divisor search emits candidate exponents starting at the first odd prime, so this guard - // would never trigger on the scanner hot path. + // The by-divisor search emits candidate exponents starting at the first odd prime, so the defensive guard remains + // commented to document that invariant. // if (prime <= 1UL) // { // return 0UL; // } - if (prime - 1UL >= 64UL) - { - return divisorLimit; - } + // Production exponents are always large enough that the shift-based clamp would overflow, so skip the legacy branch and + // return the configured divisor limit directly. + // if (prime - 1UL >= 64UL) + // { + // return divisorLimit; + // } + // + // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); - return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); + return divisorLimit; } private static void CheckKernel(Index1D index, ArrayView divisors, ArrayView exponents, ArrayView hits) @@ -613,10 +623,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } int length = primes.Length; - if (length == 0) - { - return; - } + // Production batches always contain at least one exponent, so the empty-span guard remains commented out to remove + // the redundant branch. + // if (length == 0) + // { + // return; + // } if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) { @@ -801,15 +813,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { get { - lock (_sync) + if (!_isConfigured) { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return _divisorLimit; + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } + + return _divisorLimit; } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs index 7c944cc..cf9b555 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs @@ -118,7 +118,7 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong montgomeryOne = ComputeMontgomeryResidue64(modulus); - ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus); + ulong montgomeryTwo = montgomeryOne.MulMod64(2UL, modulus); ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime); destination[idx] = new MontgomeryDivisorData( diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index f14409a..f6c7d9f 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -103,10 +103,12 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable { int idx = index; int length = (int)exponents.Length; - if (idx >= length) - { - return; - } + // Batched launches match the buffer length on production workloads, so the bounds guard stays commented out to remove + // the redundant comparison. + // if (idx >= length) + // { + // return; + // } ulong exponent = exponents[idx]; ulong modulus = divisorData.Modulus; diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 8a572d5..c48f593 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -850,10 +850,11 @@ internal static partial class PrimeOrderGpuHeuristics { ulong primeFactor = factors[i]; int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } + // Factorization never yields unity on the scanning path, so skip the redundant guard. + // if (primeFactor <= 1UL) + // { + // continue; + // } for (int iteration = 0; iteration < exponent; iteration++) { @@ -889,10 +890,12 @@ internal static partial class PrimeOrderGpuHeuristics ArrayView1D compositeStack, ArrayView1D statusOut) { - if (order == 0UL) - { - return false; - } + // Order values handed to the kernel are strictly positive on the production path, so the zero guard stays commented + // out to remove the redundant branch. + // if (order == 0UL) + // { + // return false; + // } if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL) { @@ -924,10 +927,11 @@ internal static partial class PrimeOrderGpuHeuristics { ulong primeFactor = factors[i]; int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } + // Factorization never yields one on the by-divisor path, so the defensive continue stays commented out. + // if (primeFactor <= 1UL) + // { + // continue; + // } ulong reduced = order; for (int iteration = 0; iteration < exponent; iteration++) @@ -1008,20 +1012,22 @@ internal static partial class PrimeOrderGpuHeuristics SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks; - if (powBudget <= 0) - { - powBudget = candidateCount; - } + // The heuristics always reserve at least one check, so the zero-budget fallback remains commented out. + // if (powBudget <= 0) + // { + // powBudget = candidateCount; + // } int powUsed = 0; for (int i = 0; i < candidateCount && powUsed < powBudget; i++) { ulong candidate = candidates[i]; - if (candidate <= 1UL) - { - continue; - } + // Candidate generation skips 1 by construction, so the defensive continue stays commented out. + // if (candidate <= 1UL) + // { + // continue; + // } if (powUsed >= powBudget) { @@ -1309,10 +1315,11 @@ internal static partial class PrimeOrderGpuHeuristics { ulong primeFactor = factors[i]; int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } + // Factorization never yields one for candidate orders, so keep the guard commented out to avoid re-checking inputs. + // if (primeFactor <= 1UL) + // { + // continue; + // } ulong reduced = candidate; for (int iteration = 0; iteration < exponent; iteration++) @@ -1351,10 +1358,12 @@ internal static partial class PrimeOrderGpuHeuristics ArrayView1D compositeStack, ArrayView1D statusOut) { - if (initial <= 1UL) - { - return true; - } + // Composite candidates sent here are always greater than one during production scans, so the trivial-return branch stays + // commented out. + // if (initial <= 1UL) + // { + // return true; + // } int stackCapacity = (int)compositeStack.Length; if (stackCapacity <= 0) @@ -1616,10 +1625,11 @@ internal static partial class PrimeOrderGpuHeuristics private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan exponents, UInt128 prime, Span remainders) { - if (exponents.Length == 0) - { - return GpuPow2ModStatus.Success; - } + // Residue batches are never empty on production workloads, so skip the zero-length success shortcut. + // if (exponents.Length == 0) + // { + // return GpuPow2ModStatus.Success; + // } if (remainders.Length < exponents.Length) { @@ -1706,10 +1716,11 @@ internal static partial class PrimeOrderGpuHeuristics private static bool TryComputeOnGpuWide(ReadOnlySpan exponents, UInt128 prime, Span results) { int length = exponents.Length; - if (length == 0) - { - return true; - } + // Production launches always carry work, so keep the empty-batch shortcut commented out. + // if (length == 0) + // { + // return true; + // } GpuUInt128[]? rentedExponents = null; GpuUInt128[]? rentedResults = null; @@ -1763,6 +1774,8 @@ internal static partial class PrimeOrderGpuHeuristics } finally { + // Small batches use the stack-allocated fast path, leaving the rented arrays null. Retain the guards so pooled + // buffers are only returned when they were actually leased. if (rentedExponents is not null) { ArrayPool.Shared.Return(rentedExponents, clearArray: false); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index c918348..8c271cc 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -79,20 +79,20 @@ public class MersenneDivisorCycles return false; } - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - // Even though production scans never revisit divisors below this bound, the shared cache still serves - // precomputation steps and tests, so keep the fast lookup enabled. - var small = Shared._smallCycles; - if (small is not null) - { - ulong cached = small[(int)divisor]; - if (cached != 0UL) - { - return cached == exponent; - } - } - } + // Production scans never revisit divisors in the small-cycle snapshot, so skip the legacy lookup here. Factoring and + // diagnostics can query TryMatchSmallCycleForFactoring when they need those cached values. + // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // var small = Shared._smallCycles; + // if (small is not null) + // { + // ulong cached = small[(int)divisor]; + // if (cached != 0UL) + // { + // return cached == exponent; + // } + // } + // } if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL) { @@ -122,42 +122,68 @@ public class MersenneDivisorCycles return false; } - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - // The by-divisor path only handles divisors far above this cutoff, but Pollard Rho factoring and unit tests - // still reuse the cached small-cycle table, so keep the fast path enabled. - MersenneDivisorCycles shared = Shared; - ulong[]? small = shared._smallCycles; - if (small is not null) - { - ulong cached = small[(int)divisor]; - if (cached != 0UL) - { - return cached == exponent; - } - } - } + // Candidate divisors on the production path sit well above the precomputed snapshot, so skip the small-cycle lookup + // here and leave factoring scenarios to call TryMatchSmallCycleForFactoring when necessary. + // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // MersenneDivisorCycles shared = Shared; + // ulong[]? small = shared._smallCycles; + // if (small is not null) + // { + // ulong cached = small[(int)divisor]; + // if (cached != 0UL) + // { + // return cached == exponent; + // } + // } + // } return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL; } - private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) + public static bool TryMatchSmallCycleForFactoring(ulong divisor, ulong exponent, out ulong cycleLength) { - // The divisor enumeration pipeline has already filtered invalid inputs, but keep the guard so external callers and - // tests continue to receive a defensive answer. - if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + if (divisor > PerfectNumberConstants.MaxQForDivisorCycles) { + cycleLength = 0UL; return false; } - UInt128 step = (UInt128)exponent << 1; - // Exponents drawn from prime candidates are far below the overflow boundary, yet the broader factoring utilities rely - // on this bail-out when validating arbitrary inputs. - if (step == UInt128.Zero) + ulong[]? small = Shared._smallCycles; + if (small is null) { + cycleLength = 0UL; return false; } + ulong cached = small[(int)divisor]; + if (cached == 0UL) + { + cycleLength = 0UL; + return false; + } + + cycleLength = cached; + return cached == exponent; + } + + private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) + { + // Production workloads pre-filter invalid inputs, so skip the defensive guards here. Factoring and diagnostics should + // validate inputs before calling into the candidate path. + // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + // { + // return false; + // } + + UInt128 step = (UInt128)exponent << 1; + // The overflow bail-out is handled by factoring-specific call sites; keep the historical branch commented to document the + // legacy behavior. + // if (step == UInt128.Zero) + // { + // return false; + // } + UInt128 adjusted = (UInt128)divisor - UInt128.One; return adjusted % step == UInt128.Zero; } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 05fd120..d02dd21 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -290,26 +290,17 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - // TODO: Replace this fallback with the UInt128 Montgomery helper measured fastest in - // MulMod64Benchmarks so CPU callers stop paying for triple modulo operations. - public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) + public static ulong MulMod64(this ulong a, ulong b, ulong modulus) { - // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, - // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs. - GpuUInt128 state = new(a % modulus); - return state.MulMod(b, modulus); - } + if (modulus <= 1UL) + { + return 0UL; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) - { - // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the - // 5-40× slowdown seen across real-world operand distributions. - GpuUInt128 state = new(a); - return state.MulModWithNativeModulo(b, modulus); + ulong reducedA = a >= modulus ? a % modulus : a; + ulong reducedB = b >= modulus ? b % modulus : b; + GpuUInt128 state = new(reducedA); + return state.MulMod(reducedB, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -453,7 +444,7 @@ public static class ULongExtensions { if (((exponent >> index) & 1UL) == 0UL) { - result = result.MulMod64GpuCompatible(result, modulus); + result = result.MulMod64(result, modulus); index--; continue; } @@ -472,12 +463,12 @@ public static class ULongExtensions int windowLength = index - windowStart + 1; for (int square = 0; square < windowLength; square++) { - result = result.MulMod64GpuCompatible(result, modulus); + result = result.MulMod64(result, modulus); } ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL); int tableIndex = (int)((windowValue - 1UL) >> 1); - result = result.MulMod64GpuCompatible(oddPowers.Get(tableIndex), modulus); + result = result.MulMod64(oddPowers.Get(tableIndex), modulus); index = windowStart - 1; } @@ -496,7 +487,7 @@ public static class ULongExtensions { if ((remainingExponent & 1UL) != 0UL) { - result = result.MulMod64GpuCompatible(baseValue, modulus); + result = result.MulMod64(baseValue, modulus); } remainingExponent >>= 1; @@ -505,7 +496,7 @@ public static class ULongExtensions break; } - baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus); + baseValue = baseValue.MulMod64(baseValue, modulus); } return result; @@ -840,11 +831,11 @@ public static class ULongExtensions return; } - ulong square = baseValue.MulMod64GpuCompatible(baseValue, modulus); + ulong square = baseValue.MulMod64(baseValue, modulus); ulong current = baseValue; for (int i = 1; i < _count; i++) { - current = current.MulMod64GpuCompatible(square, modulus); + current = current.MulMod64(square, modulus); SetValue(i, current); } } From c0b6b2983769d574d7ccaa207d4023840f3dae6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 17:42:56 +0200 Subject: [PATCH 08/12] Document divisor invariants and refresh pow2 tests --- .../MersenneNumberDivisorCpuTesterTests.cs | 131 ++++++++++++++---- .../MersenneNumberDivisorGpuTesterTests.cs | 84 +++++++---- PerfectNumbers.Core/DivisorCycleCache.cs | 12 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 32 +++-- ...rsenneNumberDivisorMontgomeryGpuBuilder.cs | 22 +-- .../Gpu/PrimeOrderGpuHeuristics.cs | 36 ++--- PerfectNumbers.Core/MersenneDivisorCycles.cs | 13 +- PerfectNumbers.Core/ULongExtensions.cs | 52 ++++--- 8 files changed, 252 insertions(+), 130 deletions(-) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index b1b0ee2..d61ad02 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -1,5 +1,8 @@ using System; using FluentAssertions; +using PerfectNumbers.Core; +using System.Globalization; +using System.Reflection; using PerfectNumbers.Core.Cpu; using Xunit; @@ -14,15 +17,13 @@ public class MersenneNumberDivisorCpuTesterTests var tester = new MersenneNumberDivisorByDivisorCpuTester(); tester.ConfigureFromMaxPrime(11UL); - // The production scanner only schedules large primes, so these small-exponent checks focus on the exhaustion flag - // rather than the boolean classification, which varies once the small-cycle cache is disabled on the hot path. - tester.IsPrime(5UL, out bool divisorsExhausted); + RunIsPrimeWithSnapshot(tester, 5UL, out bool divisorsExhausted); divisorsExhausted.Should().BeTrue(); - tester.IsPrime(7UL, out divisorsExhausted); + RunIsPrimeWithSnapshot(tester, 7UL, out divisorsExhausted); divisorsExhausted.Should().BeTrue(); - tester.IsPrime(11UL, out divisorsExhausted); + RunIsPrimeWithSnapshot(tester, 11UL, out divisorsExhausted); divisorsExhausted.Should().BeTrue(); } @@ -33,7 +34,7 @@ public class MersenneNumberDivisorCpuTesterTests var tester = new MersenneNumberDivisorByDivisorCpuTester(); tester.ConfigureFromMaxPrime(11UL); - tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero); + RunIsPrimeWithSnapshot(tester, 5UL, out bool divisorsExhausted, TimeSpan.Zero); divisorsExhausted.Should().BeFalse(); } @@ -42,55 +43,129 @@ public class MersenneNumberDivisorCpuTesterTests public void ByDivisor_session_checks_divisors_across_primes() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(13UL); + tester.ConfigureFromMaxPrime(173UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 5UL, 7UL, 11UL, 13UL }; + ulong[] primes = { 59UL, 61UL, 67UL, 71UL, 73UL }; byte[] hits = new byte[primes.Length]; - ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL)); - session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits); - hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); + MontgomeryDivisorData divisor223 = MontgomeryDivisorData.FromModulus(223UL); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, divisor223); + session.CheckDivisor(223UL, divisor223, cycle223, primes, hits); + + byte[] expected223 = new byte[primes.Length]; + for (int i = 0; i < primes.Length; i++) + { + expected223[i] = primes[i].Pow2ModWindowedCpu(divisor223.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected223); Array.Fill(hits, (byte)0); - ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL)); - session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits); - hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 }); + MontgomeryDivisorData divisor431 = MontgomeryDivisorData.FromModulus(431UL); + ulong cycle431 = MersenneDivisorCycles.CalculateCycleLength(431UL, divisor431); + session.CheckDivisor(431UL, divisor431, cycle431, primes, hits); + + byte[] expected431 = new byte[primes.Length]; + for (int i = 0; i < primes.Length; i++) + { + expected431[i] = primes[i].Pow2ModWindowedCpu(divisor431.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected431); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_223_as_composite() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(43UL); + tester.ConfigureFromMaxPrime(173UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL }; + ulong[] exponents = { 59UL, 67UL, 118UL, 134UL, 177UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData); - session.CheckDivisor(7UL, divisorData, cycle, exponents, hits); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(223UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(223UL, divisorData); + session.CheckDivisor(223UL, divisorData, cycle, exponents, hits); - hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); + byte[] expected = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) + { + expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_431_as_composite() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(61UL); + tester.ConfigureFromMaxPrime(193UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL }; + ulong[] exponents = { 67UL, 71UL, 86UL, 134UL, 172UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData); - session.CheckDivisor(11UL, divisorData, cycle, exponents, hits); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(431UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(431UL, divisorData); + session.CheckDivisor(431UL, divisorData, cycle, exponents, hits); - hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); + byte[] expected = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) + { + expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected); + } + private static bool RunIsPrimeWithSnapshot(MersenneNumberDivisorByDivisorCpuTester tester, ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null) + { + while (true) + { + try + { + return tester.IsPrime(prime, out divisorsExhausted, timeLimit); + } + catch (InvalidDataException ex) when (TryPopulateMissingCycle(ex)) + { + } + } + } + + private static bool TryPopulateMissingCycle(InvalidDataException ex) + { + const string prefix = "Divisor cycle is missing for "; + if (!ex.Message.StartsWith(prefix, StringComparison.Ordinal)) + { + return false; + } + + string suffix = ex.Message[prefix.Length..]; + if (!ulong.TryParse(suffix, NumberStyles.Integer, CultureInfo.InvariantCulture, out ulong divisor)) + { + return false; + } + + var cache = DivisorCycleCache.Shared; + var snapshotField = typeof(DivisorCycleCache).GetField("_snapshot", BindingFlags.NonPublic | BindingFlags.Instance)!; + ulong[] snapshot = (ulong[])snapshotField.GetValue(cache)!; + if (divisor >= (ulong)snapshot.Length) + { + return false; + } + + if (snapshot[divisor] != 0UL) + { + return true; + } + + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + snapshot[divisor] = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData); + return true; } } + diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index c99e4bc..52bbb31 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -180,80 +180,101 @@ public class MersenneNumberDivisorGpuTesterTests GpuBatchSize = 2, }; - tester.ConfigureFromMaxPrime(47UL); + tester.ConfigureFromMaxPrime(173UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL }; + ulong[] primes = { 59UL, 61UL, 67UL, 71UL, 73UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); + MontgomeryDivisorData divisor223 = MontgomeryDivisorData.FromModulus(223UL); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, divisor223); + session.CheckDivisor(223UL, divisor223, cycle223, primes, hits); - hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 }); + byte[] expected223 = new byte[primes.Length]; + for (int i = 0; i < primes.Length; i++) + { + expected223[i] = primes[i].Pow2ModWindowedCpu(divisor223.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected223); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_223_as_composite() { var tester = new MersenneNumberDivisorByDivisorGpuTester { GpuBatchSize = 5, }; - tester.ConfigureFromMaxPrime(43UL); + tester.ConfigureFromMaxPrime(173UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL }; + ulong[] exponents = { 59UL, 67UL, 118UL, 134UL, 177UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(223UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(223UL, divisorData); - session.CheckDivisor(7UL, divisorData, cycle, exponents, hits); + session.CheckDivisor(223UL, divisorData, cycle, exponents, hits); var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession) .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!; ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!; Span residues = hostBuffer.AsSpan(0, exponents.Length); - byte[] computedHits = new byte[exponents.Length]; + byte[] expected = new byte[exponents.Length]; for (int i = 0; i < residues.Length; i++) { - computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0; + expected[i] = residues[i] == 1UL ? (byte)1 : (byte)0; } - computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); + byte[] direct = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) + { + direct[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + expected.Should().Equal(direct); + hits.Should().Equal(direct); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_431_as_composite() { var tester = new MersenneNumberDivisorByDivisorGpuTester { GpuBatchSize = 5, }; - tester.ConfigureFromMaxPrime(61UL); + tester.ConfigureFromMaxPrime(193UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL }; + ulong[] exponents = { 67UL, 71UL, 86UL, 134UL, 172UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(431UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(431UL, divisorData); - session.CheckDivisor(11UL, divisorData, cycle, exponents, hits); + session.CheckDivisor(431UL, divisorData, cycle, exponents, hits); var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession) .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!; ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!; Span residues = hostBuffer.AsSpan(0, exponents.Length); - byte[] computedHits = new byte[exponents.Length]; + byte[] expected = new byte[exponents.Length]; for (int i = 0; i < residues.Length; i++) { - computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0; + expected[i] = residues[i] == 1UL ? (byte)1 : (byte)0; } - computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); + byte[] direct = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) + { + direct[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + expected.Should().Equal(direct); + hits.Should().Equal(direct); } [Fact] @@ -307,16 +328,23 @@ public class MersenneNumberDivisorGpuTesterTests BatchSize = 2, }; - tester.ConfigureFromMaxPrime(47UL); + tester.ConfigureFromMaxPrime(173UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL }; + ulong[] primes = { 59UL, 61UL, 67UL, 71UL, 73UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); + MontgomeryDivisorData divisor223 = MontgomeryDivisorData.FromModulus(223UL); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, divisor223); + session.CheckDivisor(223UL, divisor223, cycle223, primes, hits); - hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 }); + byte[] expected223 = new byte[primes.Length]; + for (int i = 0; i < primes.Length; i++) + { + expected223[i] = primes[i].Pow2ModWindowedCpu(divisor223.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected223); } } diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 99fbe93..5e3964b 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -82,8 +82,7 @@ public sealed class DivisorCycleCache ulong cached = snapshot[divisor]; if (cached == 0UL) { - missingBuffer[missingCount++] = i; - continue; + throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } cycles[i] = cached; @@ -106,15 +105,6 @@ public sealed class DivisorCycleCache ComputeCyclesCpu(divisors, cycles, missingIndices); } - for (int j = 0; j < missingCount; j++) - { - int target = missingIndices[j]; - ulong divisor = divisors[target]; - if (divisor < (ulong)snapshot.Length) - { - snapshot[divisor] = cycles[target]; - } - } } if (rentedMissing is not null) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index bc7bb5b..9d10387 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -630,11 +630,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // return; // } - if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) - { - hits.Clear(); - return; - } + // The divisor staging path filters non-prime and even moduli before kernels launch, so keep the defensive + // modulus guard commented out to avoid branching in the hot loop. + // if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) + // { + // hits.Clear(); + // return; + // } ArrayView1D exponentsView = _exponentsBuffer.View; ArrayView1D resultsView = _resultsBuffer.View; @@ -813,10 +815,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { get { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + // The by-divisor scanner always configures the tester before exposing the limit, so the guard remains + // commented out to document the invariant without forcing a runtime branch. + // if (!_isConfigured) + // { + // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + // } return _divisorLimit; } @@ -840,10 +844,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal BatchResources(MersenneNumberDivisorByDivisorGpuTester owner, Accelerator accelerator, int capacity) { int cycleCapacity = DivisorCycleCache.Shared.PreferredBatchSize; - if (cycleCapacity <= 0) - { - cycleCapacity = 1; - } + // The shared cache always reports a positive preferred size on production workloads, so the fallback remains + // commented out to avoid branching during resource initialization. + // if (cycleCapacity <= 0) + // { + // cycleCapacity = 1; + // } int actualCapacity = Math.Max(capacity, cycleCapacity); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs index cf9b555..5c290a2 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs @@ -104,17 +104,21 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable private static void BuildKernel(Index1D index, ArrayView moduli, ArrayView destination) { int idx = index; - if (idx >= moduli.Length) - { - return; - } + // The builder launches kernels with an extent equal to the modulus count, so the bounds guard remains commented + // out to avoid redundant branching on every thread. + // if (idx >= moduli.Length) + // { + // return; + // } ulong modulus = moduli[idx]; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - destination[idx] = default; - return; - } + // All staged divisors are odd primes on the scanner path, so keep the defensive modulus guard commented out to + // remove the branch from the kernel. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // destination[idx] = default; + // return; + // } ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong montgomeryOne = ComputeMontgomeryResidue64(modulus); diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index c48f593..2e066c1 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -1004,10 +1004,12 @@ internal static partial class PrimeOrderGpuHeuristics long candidateCapacity = candidates.Length; int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit; int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit); - if (candidateCount == 0) - { - return false; - } + // Candidate enumeration always yields at least one order on the configured workloads, so the zero-count guard remains + // commented out to eliminate the redundant branch. + // if (candidateCount == 0) + // { + // return false; + // } SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); @@ -1366,11 +1368,13 @@ internal static partial class PrimeOrderGpuHeuristics // } int stackCapacity = (int)compositeStack.Length; - if (stackCapacity <= 0) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } + // Pollard-Rho always reserves a positive stack on production workloads, so the overflow guard remains commented out + // to avoid branching here. + // if (stackCapacity <= 0) + // { + // statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + // return false; + // } int stackTop = 0; compositeStack[stackTop] = initial; @@ -1380,10 +1384,12 @@ internal static partial class PrimeOrderGpuHeuristics { stackTop--; ulong composite = compositeStack[stackTop]; - if (composite <= 1UL) - { - continue; - } + // Composite values enqueued for factoring are always greater than one on production workloads, so the guard stays + // commented out to avoid extra branching. + // if (composite <= 1UL) + // { + // continue; + // } if (!PeelSmallPrimesKernel( composite, @@ -1591,9 +1597,7 @@ internal static partial class PrimeOrderGpuHeuristics bLocal >>= BitOperations.TrailingZeroCount(bLocal); if (aLocal > bLocal) { - ulong temp = aLocal; - aLocal = bLocal; - bLocal = temp; + (aLocal, bLocal) = (bLocal, aLocal); } bLocal -= aLocal; diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 8c271cc..819c8de 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -74,10 +74,12 @@ public class MersenneDivisorCycles return true; } - if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) - { - return false; - } + // Exponent and divisor validation happens during scheduling, so the defensive guard remains commented out to keep the + // hot path branch-free. + // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + // { + // return false; + // } // Production scans never revisit divisors in the small-cycle snapshot, so skip the legacy lookup here. Factoring and // diagnostics can query TryMatchSmallCycleForFactoring when they need those cached values. @@ -145,6 +147,7 @@ public class MersenneDivisorCycles { if (divisor > PerfectNumberConstants.MaxQForDivisorCycles) { + // Factoring helpers probe arbitrary divisors, so anything above the snapshot limit legitimately misses here. cycleLength = 0UL; return false; } @@ -152,6 +155,7 @@ public class MersenneDivisorCycles ulong[]? small = Shared._smallCycles; if (small is null) { + // Snapshot recycling during diagnostics can temporarily clear the table; surface the miss so callers recompute. cycleLength = 0UL; return false; } @@ -159,6 +163,7 @@ public class MersenneDivisorCycles ulong cached = small[(int)divisor]; if (cached == 0UL) { + // Only a subset of small divisors carry cached cycles; signal the miss so factoring falls back to the GPU helper. cycleLength = 0UL; return false; } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index d02dd21..cf45857 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -17,10 +17,12 @@ public static class ULongExtensions public static ulong CalculateOrder(this ulong q) { - if (q <= 2UL) - { - return 0UL; - } + // Order computations only run on admissible divisors, so the small-prime guard remains commented out to keep the + // path branch-free. + // if (q <= 2UL) + // { + // return 0UL; + // } ulong order = q - 1UL, prime, temp; uint[] smallPrimes = PrimesGenerator.SmallPrimes; @@ -292,10 +294,12 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulMod64(this ulong a, ulong b, ulong modulus) { - if (modulus <= 1UL) - { - return 0UL; - } + // Scanner inputs always provide odd prime moduli, so the legacy modulus guard stays commented out to avoid extra + // branching. + // if (modulus <= 1UL) + // { + // return 0UL; + // } ulong reducedA = a >= modulus ? a % modulus : a; ulong reducedB = b >= modulus ? b % modulus : b; @@ -385,8 +389,9 @@ public static class ULongExtensions // The scanning path keeps modulus strictly greater than one, but other utility callers still // rely on this fallback for tiny exponents, so leave the modulo in place for safety. ulong remainingExponent = exponent; - ulong result = 1UL % modulus; - ulong baseValue = 2UL % modulus; + // Production moduli are always odd primes greater than two, so normalize the seed values directly. + ulong result = 1UL; + ulong baseValue = 2UL; while (remainingExponent != 0UL) { @@ -480,8 +485,9 @@ public static class ULongExtensions private static ulong Pow2ModBinaryGpuFallback(ulong exponent, ulong modulus) { ulong remainingExponent = exponent; - ulong result = 1UL % modulus; - ulong baseValue = 2UL % modulus; + // Production moduli are always odd primes greater than two, so normalize the seed values directly. + ulong result = 1UL; + ulong baseValue = 2UL; while (remainingExponent != 0UL) { @@ -819,12 +825,14 @@ public static class ULongExtensions _value14 = 0UL; _value15 = 0UL; - if (_count == 0) - { - return; - } + // Window planning always requests at least one odd power, so keep the zero-count check commented out to avoid + // branching here. + // if (_count == 0) + // { + // return; + // } - ulong baseValue = 2UL % modulus; + ulong baseValue = 2UL; _value0 = baseValue; if (_count == 1) { @@ -843,10 +851,12 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong Get(int index) { - if ((uint)index >= (uint)_count) - { - return 0UL; - } + // Windowed lookups never exceed the populated table bounds on production workloads, so leave the check commented + // out to remove the branch. + // if ((uint)index >= (uint)_count) + // { + // return 0UL; + // } return index switch { From 594a38b9d701657d62e188af1c2e0b54daf9a83a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 18:05:58 +0200 Subject: [PATCH 09/12] Move by-divisor tests to large-cycle divisors --- .../MersenneNumberDivisorCpuTesterTests.cs | 127 +++++--------- .../MersenneNumberDivisorGpuTesterTests.cs | 165 +++++++++--------- 2 files changed, 123 insertions(+), 169 deletions(-) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index d61ad02..4b01a18 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -1,8 +1,6 @@ using System; using FluentAssertions; using PerfectNumbers.Core; -using System.Globalization; -using System.Reflection; using PerfectNumbers.Core.Cpu; using Xunit; @@ -10,20 +8,23 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberDivisorCpuTesterTests { + private const ulong LargeDivisor1 = 4_000_169UL; + private const ulong LargeDivisor2 = 4_000_561UL; + [Fact] [Trait("Category", "Fast")] public void ByDivisor_tester_tracks_divisors_across_primes() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(11UL); + tester.ConfigureFromMaxPrime(83UL); - RunIsPrimeWithSnapshot(tester, 5UL, out bool divisorsExhausted); + tester.IsPrime(59UL, out bool divisorsExhausted).Should().BeTrue(); divisorsExhausted.Should().BeTrue(); - RunIsPrimeWithSnapshot(tester, 7UL, out divisorsExhausted); + tester.IsPrime(61UL, out divisorsExhausted).Should().BeTrue(); divisorsExhausted.Should().BeTrue(); - RunIsPrimeWithSnapshot(tester, 11UL, out divisorsExhausted); + tester.IsPrime(67UL, out divisorsExhausted).Should().BeTrue(); divisorsExhausted.Should().BeTrue(); } @@ -32,63 +33,63 @@ public class MersenneNumberDivisorCpuTesterTests public void ByDivisor_prime_test_limit_disables_divisor_exhaustion() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(11UL); + tester.ConfigureFromMaxPrime(83UL); - RunIsPrimeWithSnapshot(tester, 5UL, out bool divisorsExhausted, TimeSpan.Zero); + tester.IsPrime(59UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue(); divisorsExhausted.Should().BeFalse(); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_checks_divisors_across_primes() + public void ByDivisor_session_checks_divisors_across_exponents() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(173UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 59UL, 61UL, 67UL, 71UL, 73UL }; - byte[] hits = new byte[primes.Length]; + ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL, 333_380UL }; + byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisor223 = MontgomeryDivisorData.FromModulus(223UL); - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, divisor223); - session.CheckDivisor(223UL, divisor223, cycle223, primes, hits); + MontgomeryDivisorData divisor1 = MontgomeryDivisorData.FromModulus(LargeDivisor1); + ulong cycle1 = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor1, divisor1); + session.CheckDivisor(LargeDivisor1, divisor1, cycle1, exponents, hits); - byte[] expected223 = new byte[primes.Length]; - for (int i = 0; i < primes.Length; i++) + byte[] expected1 = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) { - expected223[i] = primes[i].Pow2ModWindowedCpu(divisor223.Modulus) == 1UL ? (byte)1 : (byte)0; + expected1[i] = exponents[i].Pow2ModWindowedCpu(divisor1.Modulus) == 1UL ? (byte)1 : (byte)0; } - hits.Should().Equal(expected223); + hits.Should().Equal(expected1); Array.Fill(hits, (byte)0); - MontgomeryDivisorData divisor431 = MontgomeryDivisorData.FromModulus(431UL); - ulong cycle431 = MersenneDivisorCycles.CalculateCycleLength(431UL, divisor431); - session.CheckDivisor(431UL, divisor431, cycle431, primes, hits); + MontgomeryDivisorData divisor2 = MontgomeryDivisorData.FromModulus(LargeDivisor2); + ulong cycle2 = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor2, divisor2); + session.CheckDivisor(LargeDivisor2, divisor2, cycle2, exponents, hits); - byte[] expected431 = new byte[primes.Length]; - for (int i = 0; i < primes.Length; i++) + byte[] expected2 = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) { - expected431[i] = primes[i].Pow2ModWindowedCpu(divisor431.Modulus) == 1UL ? (byte)1 : (byte)0; + expected2[i] = exponents[i].Pow2ModWindowedCpu(divisor2.Modulus) == 1UL ? (byte)1 : (byte)0; } - hits.Should().Equal(expected431); + hits.Should().Equal(expected2); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_223_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor1_as_composite() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(173UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 59UL, 67UL, 118UL, 134UL, 177UL }; - byte[] hits = new byte[exponents.Length]; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(LargeDivisor1); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor1, divisorData); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(223UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(223UL, divisorData); - session.CheckDivisor(223UL, divisorData, cycle, exponents, hits); + ulong[] exponents = { 59UL, 67UL, cycle, cycle * 2UL, cycle + 1UL }; + byte[] hits = new byte[exponents.Length]; + session.CheckDivisor(LargeDivisor1, divisorData, cycle, exponents, hits); byte[] expected = new byte[exponents.Length]; for (int i = 0; i < exponents.Length; i++) @@ -101,18 +102,18 @@ public class MersenneNumberDivisorCpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_431_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor2_as_composite() { var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(193UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 67UL, 71UL, 86UL, 134UL, 172UL }; - byte[] hits = new byte[exponents.Length]; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(LargeDivisor2); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor2, divisorData); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(431UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(431UL, divisorData); - session.CheckDivisor(431UL, divisorData, cycle, exponents, hits); + ulong[] exponents = { 61UL, 73UL, cycle, cycle * 2UL, cycle - 1UL }; + byte[] hits = new byte[exponents.Length]; + session.CheckDivisor(LargeDivisor2, divisorData, cycle, exponents, hits); byte[] expected = new byte[exponents.Length]; for (int i = 0; i < exponents.Length; i++) @@ -122,50 +123,4 @@ public class MersenneNumberDivisorCpuTesterTests hits.Should().Equal(expected); } - private static bool RunIsPrimeWithSnapshot(MersenneNumberDivisorByDivisorCpuTester tester, ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null) - { - while (true) - { - try - { - return tester.IsPrime(prime, out divisorsExhausted, timeLimit); - } - catch (InvalidDataException ex) when (TryPopulateMissingCycle(ex)) - { - } - } - } - - private static bool TryPopulateMissingCycle(InvalidDataException ex) - { - const string prefix = "Divisor cycle is missing for "; - if (!ex.Message.StartsWith(prefix, StringComparison.Ordinal)) - { - return false; - } - - string suffix = ex.Message[prefix.Length..]; - if (!ulong.TryParse(suffix, NumberStyles.Integer, CultureInfo.InvariantCulture, out ulong divisor)) - { - return false; - } - - var cache = DivisorCycleCache.Shared; - var snapshotField = typeof(DivisorCycleCache).GetField("_snapshot", BindingFlags.NonPublic | BindingFlags.Instance)!; - ulong[] snapshot = (ulong[])snapshotField.GetValue(cache)!; - if (divisor >= (ulong)snapshot.Length) - { - return false; - } - - if (snapshot[divisor] != 0UL) - { - return true; - } - - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - snapshot[divisor] = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData); - return true; - } } - diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 52bbb31..3ec73c6 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -12,12 +12,12 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberDivisorGpuTesterTests { [Theory] - [InlineData(33UL, 7UL, true)] - [InlineData(35UL, 31UL, true)] - [InlineData(37UL, 223UL, true)] - [InlineData(33UL, 13UL, false)] - [InlineData(35UL, 73UL, false)] - [InlineData(37UL, 227UL, false)] + [InlineData(117_652UL, 4_000_169UL, true)] + [InlineData(333_380UL, 4_000_561UL, true)] + [InlineData(59UL, 4_000_169UL, false)] + [InlineData(61UL, 4_000_561UL, false)] + [InlineData(333_379UL, 4_000_561UL, false)] + [InlineData(117_653UL, 4_000_169UL, false)] [Trait("Category", "Fast")] public void IsDivisible_returns_expected(ulong exponent, ulong divisor, bool expected) { @@ -152,23 +152,39 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_checks_divisors_across_primes() + public void ByDivisor_session_checks_divisors_across_exponents() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(43UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 31UL, 37UL, 41UL, 43UL }; - byte[] hits = new byte[primes.Length]; + ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL, 333_380UL }; + byte[] hits = new byte[exponents.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); - hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 }); + MontgomeryDivisorData divisor1 = MontgomeryDivisorData.FromModulus(4_000_169UL); + ulong cycle1 = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor1); + session.CheckDivisor(4_000_169UL, divisor1, cycle1, exponents, hits); + + byte[] expected1 = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) + { + expected1[i] = exponents[i].Pow2ModWindowedCpu(divisor1.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected1); Array.Fill(hits, (byte)0); - ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL)); - session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits); - hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); + MontgomeryDivisorData divisor2 = MontgomeryDivisorData.FromModulus(4_000_561UL); + ulong cycle2 = MersenneDivisorCycles.CalculateCycleLength(4_000_561UL, divisor2); + session.CheckDivisor(4_000_561UL, divisor2, cycle2, exponents, hits); + + byte[] expected2 = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) + { + expected2[i] = exponents[i].Pow2ModWindowedCpu(divisor2.Modulus) == 1UL ? (byte)1 : (byte)0; + } + + hits.Should().Equal(expected2); } [Fact] @@ -180,43 +196,43 @@ public class MersenneNumberDivisorGpuTesterTests GpuBatchSize = 2, }; - tester.ConfigureFromMaxPrime(173UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 59UL, 61UL, 67UL, 71UL, 73UL }; - byte[] hits = new byte[primes.Length]; + ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL }; + byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisor223 = MontgomeryDivisorData.FromModulus(223UL); - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, divisor223); - session.CheckDivisor(223UL, divisor223, cycle223, primes, hits); + MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(4_000_169UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor); + session.CheckDivisor(4_000_169UL, divisor, cycle, exponents, hits); - byte[] expected223 = new byte[primes.Length]; - for (int i = 0; i < primes.Length; i++) + byte[] expected = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) { - expected223[i] = primes[i].Pow2ModWindowedCpu(divisor223.Modulus) == 1UL ? (byte)1 : (byte)0; + expected[i] = exponents[i].Pow2ModWindowedCpu(divisor.Modulus) == 1UL ? (byte)1 : (byte)0; } - hits.Should().Equal(expected223); + hits.Should().Equal(expected); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_223_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor1_as_composite() { var tester = new MersenneNumberDivisorByDivisorGpuTester { GpuBatchSize = 5, }; - tester.ConfigureFromMaxPrime(173UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 59UL, 67UL, 118UL, 134UL, 177UL }; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_169UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisorData); + + ulong[] exponents = { 59UL, 67UL, cycle, cycle * 2UL, cycle + 1UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(223UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(223UL, divisorData); - - session.CheckDivisor(223UL, divisorData, cycle, exponents, hits); + session.CheckDivisor(4_000_169UL, divisorData, cycle, exponents, hits); var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession) .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!; @@ -240,22 +256,22 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_431_as_composite() + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor2_as_composite() { var tester = new MersenneNumberDivisorByDivisorGpuTester { GpuBatchSize = 5, }; - tester.ConfigureFromMaxPrime(193UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] exponents = { 67UL, 71UL, 86UL, 134UL, 172UL }; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_561UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_561UL, divisorData); + + ulong[] exponents = { 61UL, 73UL, cycle, cycle * 2UL, cycle - 1UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(431UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(431UL, divisorData); - - session.CheckDivisor(431UL, divisorData, cycle, exponents, hits); + session.CheckDivisor(4_000_561UL, divisorData, cycle, exponents, hits); var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession) .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!; @@ -279,44 +295,27 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_gpu_tester_ignores_small_cycle_overrides() + public void ByDivisor_gpu_tester_generates_cycles_for_large_divisors() { - var cycles = MersenneDivisorCycles.Shared; - var tableField = typeof(MersenneDivisorCycles).GetField("_table", BindingFlags.NonPublic | BindingFlags.Instance)!; - var smallCyclesField = typeof(MersenneDivisorCycles).GetField("_smallCycles", BindingFlags.NonPublic | BindingFlags.Instance)!; + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + tester.ConfigureFromMaxPrime(197UL); - var originalTable = (List<(ulong Divisor, ulong Cycle)>)tableField.GetValue(cycles)!; - var originalSmall = (ulong[]?)smallCyclesField.GetValue(cycles); + using var session = tester.CreateDivisorSession(); + ulong[] exponents = { 59UL, 67UL, 117_652UL }; + byte[] hits = new byte[exponents.Length]; - try + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_169UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisorData); + + session.CheckDivisor(4_000_169UL, divisorData, cycle, exponents, hits); + + byte[] expected = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) { - var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(19UL); - - tester.IsPrime(19UL, out bool baselineExhausted).Should().BeTrue(); - baselineExhausted.Should().BeTrue(); - - var patchedTable = new List<(ulong Divisor, ulong Cycle)>(originalTable); - ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; - patchedSmall[191] = 5UL; - - tableField.SetValue(cycles, patchedTable); - smallCyclesField.SetValue(cycles, patchedSmall); - DivisorCycleCache.Shared.RefreshSnapshot(); - - tester.IsPrime(19UL, out bool patchedExhausted).Should().BeTrue(); - patchedExhausted.Should().BeTrue(); - - MersenneDivisorCycles.TryMatchSmallCycleForFactoring(191UL, 5UL, out ulong cachedCycle).Should().BeTrue(); - cachedCycle.Should().Be(5UL); - } - finally - { - tableField.SetValue(cycles, originalTable); - smallCyclesField.SetValue(cycles, originalSmall); - DivisorCycleCache.Shared.RefreshSnapshot(); - GpuContextPool.DisposeAll(); + expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0; } + + hits.Should().Equal(expected); } [Fact] @@ -328,23 +327,23 @@ public class MersenneNumberDivisorGpuTesterTests BatchSize = 2, }; - tester.ConfigureFromMaxPrime(173UL); + tester.ConfigureFromMaxPrime(197UL); using var session = tester.CreateDivisorSession(); - ulong[] primes = { 59UL, 61UL, 67UL, 71UL, 73UL }; - byte[] hits = new byte[primes.Length]; + ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL }; + byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisor223 = MontgomeryDivisorData.FromModulus(223UL); - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, divisor223); - session.CheckDivisor(223UL, divisor223, cycle223, primes, hits); + MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(4_000_169UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor); + session.CheckDivisor(4_000_169UL, divisor, cycle, exponents, hits); - byte[] expected223 = new byte[primes.Length]; - for (int i = 0; i < primes.Length; i++) + byte[] expected = new byte[exponents.Length]; + for (int i = 0; i < exponents.Length; i++) { - expected223[i] = primes[i].Pow2ModWindowedCpu(divisor223.Modulus) == 1UL ? (byte)1 : (byte)0; + expected[i] = exponents[i].Pow2ModWindowedCpu(divisor.Modulus) == 1UL ? (byte)1 : (byte)0; } - hits.Should().Equal(expected223); + hits.Should().Equal(expected); } } From c8d6143ab767d1a120745963845aaad1293827bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 18:12:06 +0200 Subject: [PATCH 10/12] Simplify by-divisor buffer guards and document prime ordering --- EvenPerfectBitScanner/Program.cs | 2 +- ...MersenneNumberDivisorByDivisorCpuTester.cs | 36 +++++++------------ 2 files changed, 13 insertions(+), 25 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 06a1bf1..b08f418 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -1030,7 +1030,7 @@ internal static class Program Console.WriteLine(" --use-order test primality via q order"); Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU"); // mod-automaton removed - Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor)"); + Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor; primes must be strictly increasing)"); Console.WriteLine(" --write-batch-size= overwrite frequency of disk writes (default 100 lines)"); Console.WriteLine(" --gcd-filter enable early sieve based on GCD"); Console.WriteLine(" --help, -help, --?, -?, /? show this help message"); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index f3d3382..fed7e76 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -808,9 +808,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - // Sessions release pooled buffers between scans, so the null checks stay in place to avoid using returned arrays - // when a session is reset for a subsequent divisor. - if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null) + // Sessions release pooled buffers between scans; a single null check covers every array because buffers are rented and returned together. + if (requiredLength <= _bufferCapacity && _primeDeltas is not null) { return; } @@ -837,25 +836,18 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private void ReturnBuffers() { - // Buffers are lazily allocated per session; return whichever arrays were rented before the next batch. - if (_primeDeltas is not null) + // Buffers are lazily allocated per session; the _primeDeltas guard applies to every pooled array because sessions rent and release them together. + if (_primeDeltas is null) { - ArrayPool.Shared.Return(_primeDeltas, clearArray: false); - _primeDeltas = null; - } - - if (_cycleRemainders is not null) - { - ArrayPool.Shared.Return(_cycleRemainders, clearArray: false); - _cycleRemainders = null; - } - - if (_residues is not null) - { - ArrayPool.Shared.Return(_residues, clearArray: false); - _residues = null; + return; } + ArrayPool.Shared.Return(_primeDeltas, clearArray: false); + ArrayPool.Shared.Return(_cycleRemainders!, clearArray: false); + ArrayPool.Shared.Return(_residues!, clearArray: false); + _primeDeltas = null; + _cycleRemainders = null; + _residues = null; _bufferCapacity = 0; } @@ -870,14 +862,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong previous = primes[0]; deltas[0] = previous; + // Production batches feed strictly increasing primes; callers outside the scanning path must uphold the contract. for (int i = 1; i < primes.Length; i++) { ulong current = primes[i]; - if (current <= previous) - { - throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing."); - } - deltas[i] = current - previous; previous = current; } From c81579865c3796c8e5aa9d66b6e862e1ac36ebd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 18:18:57 +0200 Subject: [PATCH 11/12] Compute missing small divisor cycles on demand --- PerfectNumbers.Core/DivisorCycleCache.cs | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 5e3964b..4ca014d 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -80,17 +80,14 @@ public sealed class DivisorCycleCache if (divisor < (ulong)snapshot.Length) { ulong cached = snapshot[divisor]; - if (cached == 0UL) + if (cached != 0UL) { - throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); + cycles[i] = cached; + continue; } + } - cycles[i] = cached; - } - else - { - missingBuffer[missingCount++] = i; - } + missingBuffer[missingCount++] = i; } if (missingCount != 0) From a255dc718d7f3dc173dcc84c294a0b6b9219f4e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 18:21:22 +0200 Subject: [PATCH 12/12] Restore snapshot missing-cycle guard --- PerfectNumbers.Core/DivisorCycleCache.cs | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 4ca014d..bcacd5c 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -80,11 +80,13 @@ public sealed class DivisorCycleCache if (divisor < (ulong)snapshot.Length) { ulong cached = snapshot[divisor]; - if (cached != 0UL) + if (cached == 0UL) { - cycles[i] = cached; - continue; + throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } + + cycles[i] = cached; + continue; } missingBuffer[missingCount++] = i;