diff --git a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs index d618610..090c163 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs @@ -22,12 +22,28 @@ public class GpuPow2ModBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Windowed (8-bit) exponentiation that dominated every dataset: 21.5 μs for the full-width modulus, 7.87 μs for the large + /// set, 265 ns for medium, and 82.97 ns for the small modulus sample. + /// + /// + /// Observed means: FullWidthExponentFullWidthModulus 21,481.71 ns (1.00×), LargeExponentHighWordModulus 7,874.85 ns, + /// MediumExponentPrimeModulus 264.97 ns, SmallExponentSmallModulus 82.97 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 ProcessEightBitWindows() { return GpuUInt128.Pow2Mod(Input.Exponent, Input.Modulus); } + /// + /// Bit-by-bit fallback; competitive only on medium/small moduli (266 ns / 85.39 ns) while trailing badly on the widest + /// modulus (50.97 μs, 2.37× slower) and large input (8.85 μs, 1.12× slower). + /// + /// + /// Observed means: FullWidthExponentFullWidthModulus 50,967.07 ns (2.37×), LargeExponentHighWordModulus 8,851.95 ns, + /// MediumExponentPrimeModulus 265.98 ns, SmallExponentSmallModulus 85.39 ns. + /// [Benchmark] public GpuUInt128 ProcessSingleBits() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs index 5274cb9..16efefb 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs @@ -23,6 +23,13 @@ public class GpuUInt128AddBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Materializes the carry in a local before updating the struct; stayed between 0.459 ns and 0.474 ns across all operand + /// patterns, giving it the edge on carry-heavy cases. + /// + /// + /// Observed means: AllBitsSet 0.4737 ns (1.00×), CarryIntoHigh 0.4586 ns, IncrementLow 0.4653 ns, MixedOperands 0.4628 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 CarryMaterialisedInLocal() { @@ -31,6 +38,13 @@ public class GpuUInt128AddBenchmarks return value; } + /// + /// Computes the carry inline within the expression; matches the baseline on balanced additions (0.474 ns mixed) but trails + /// slightly when the carry bubbles into the high word (0.491 ns) or chains through increments (0.508 ns). + /// + /// + /// Observed means: AllBitsSet 0.4676 ns (0.98×), CarryIntoHigh 0.4907 ns, IncrementLow 0.5082 ns, MixedOperands 0.4741 ns. + /// [Benchmark] public GpuUInt128 CarryInExpression() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs index 2fa7947..312aebb 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs @@ -26,12 +26,26 @@ public class GpuUInt128BinaryGcdBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Baseline that reuses locals between iterations; measured 195–196 ns on high-entropy/high-word-only inputs, 157.8 ns on + /// low-word-heavy pairs, and 7.32 ns on tiny operands. + /// + /// + /// Observed means: HighEntropy 195.593 ns (1.00×), HighWordOnly 196.218 ns, LowWordHeavy 157.787 ns, SmallOperands 7.317 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 ReusedVariables() { return BinaryGcdWithReusedVariables(Input.Left, Input.Right); } + /// + /// Allocates a temporary struct per loop iteration; trims 4–9% off the large inputs (186–186 ns) and 18% off small operands + /// (5.99 ns), making it the better all-around choice when GC pressure is acceptable. + /// + /// + /// Observed means: HighEntropy 186.615 ns (0.95×), HighWordOnly 185.849 ns, LowWordHeavy 143.686 ns, SmallOperands 5.987 ns. + /// [Benchmark] public GpuUInt128 TemporaryStructPerIteration() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs index d24d3ab..43d8305 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs @@ -24,15 +24,31 @@ public class GpuUInt128Montgomery64Benchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Extension-based Montgomery multiply that ran in 3.64 ns across every operand mix (dense, mixed, near-modulus, tiny), + /// making it the reference CPU helper. + /// + /// + /// Observed means: DenseOperands 3.663 ns (1.00×), MixedMagnitude 3.636 ns, NearModulus 3.651 ns, TinyOperands 3.639 ns. + /// [Benchmark(Baseline = true)] public ulong ExtensionMontgomeryMultiply() { return Input.Left.MontgomeryMultiply(Input.Right, Input.Modulus, Input.NPrime); } + /// + /// Struct-based GPU-style multiply; matches device arithmetic but costs 24.9 ns regardless of the operand pattern, + /// approximately 6.8× slower than the extension helper. + /// + /// + /// Observed means: DenseOperands 24.904 ns (6.80×), MixedMagnitude 24.886 ns, NearModulus 24.893 ns, TinyOperands 24.869 ns. + /// [Benchmark] public ulong GpuStructMontgomeryMultiply() { + // TODO: Migrate remaining callers of MulModMontgomery64 to the extension helper; the struct + // emulation is ~6.8× slower but kept for GPU parity benchmarks. GpuUInt128 state = new(Input.Left); return state.MulModMontgomery64(Input.Right, Input.Modulus, Input.NPrime, Input.R2); } diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs index 0939a84..6d125ff 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs @@ -23,6 +23,13 @@ public class GpuUInt128Mul64Benchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Baseline that materializes the high product in a local; measured 1.53–1.55 ns across all operand patterns, keeping it in + /// front of the inline variant. + /// + /// + /// Observed means: AllBitsSet 1.551 ns (1.00×), HighByLow 1.532 ns, LowWordMax 1.548 ns, TinyOperands 1.529 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 HighProductMaterializedInLocal() { @@ -31,6 +38,12 @@ public class GpuUInt128Mul64Benchmarks return value; } + /// + /// Inline variant that avoids the temporary struct field update, but costs roughly 4–5% extra work with 1.59–1.63 ns means. + /// + /// + /// Observed means: AllBitsSet 1.631 ns (1.05×), HighByLow 1.594 ns, LowWordMax 1.623 ns, TinyOperands 1.602 ns. + /// [Benchmark] public GpuUInt128 HighProductInline() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs index 100d293..742bb35 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs @@ -39,6 +39,13 @@ public class GpuUInt128MulModBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// In-place multiply-modulo; great for low-word-heavy operands (6.88 ns) and tiny pairs (7.15 ns) but notably slower on + /// large moduli (193–201 ns). + /// + /// + /// Observed means: HighWordModulus 192.794 ns (1.00×), LowWordHeavy 6.883 ns, MixedMagnitude 201.468 ns, TinyOperands 7.146 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 InPlaceMulMod() { @@ -48,6 +55,13 @@ public class GpuUInt128MulModBenchmarks return value; } + /// + /// Allocates fresh temporaries per iteration; pays off on every dataset with 4.66 ns on tiny operands, 117 ns on + /// low-word-heavy inputs, 138 ns on mixed magnitude, and 184 ns on high-word moduli. + /// + /// + /// Observed means: HighWordModulus 184.322 ns (0.96×), LowWordHeavy 116.993 ns, MixedMagnitude 138.684 ns, TinyOperands 4.659 ns. + /// [Benchmark] public GpuUInt128 AllocatePerIteration() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs index a2baa51..1f385f3 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs @@ -33,12 +33,26 @@ public class GpuUInt128MulModByLimbBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Legacy implementation that allocates intermediate limbs; delivered 10.7 ns on tiny operands, 12.7 ns on high-word-dominant + /// pairs, and 16.3 ns on mixed magnitude, making it the faster choice overall. + /// + /// + /// Observed means: HighWordDominant 12.69 ns (1.00×), MixedMagnitude 16.26 ns, TinyOperands 10.74 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 LegacyAllocating() { return MulModByLimbLegacy(Input.Left, Input.Right, Input.Modulus); } + /// + /// In-place reduction variant; avoids allocations but costs 12.9–25.1 ns, trailing the legacy path by 20–55% depending on the + /// operand mix. + /// + /// + /// Observed means: HighWordDominant 17.88 ns (1.41×), MixedMagnitude 25.14 ns, TinyOperands 12.92 ns. + /// [Benchmark] public GpuUInt128 InPlaceReduction() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs index f15cab9..535d6df 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs @@ -22,6 +22,13 @@ public class GpuUInt128NativeModuloBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Performs the modulo immediately after loading the left operand; runs in 7.28 ns on tiny inputs and 38–58 ns on dense or + /// mixed workloads, keeping it well ahead of the deferred path. + /// + /// + /// Observed means: DenseOperands 38.127 ns (1.00×), MixedMagnitude 58.555 ns, NearModulus 54.691 ns, TinyOperands 7.278 ns. + /// [Benchmark(Baseline = true)] public ulong ImmediateModulo() { @@ -29,9 +36,18 @@ public class GpuUInt128NativeModuloBenchmarks return state.MulMod(Input.Right, Input.Modulus); } + /// + /// Defers the modulo until after the multiply; only helps on tiny operands (18.15 ns) but is 4–8× slower on dense or + /// near-modulus workloads at 256–299 ns. + /// + /// + /// Observed means: DenseOperands 298.926 ns (7.83×), MixedMagnitude 255.624 ns, NearModulus 293.280 ns, TinyOperands 18.149 ns. + /// [Benchmark] public ulong DeferredNativeModulo() { + // TODO: Drop MulModWithNativeModulo from production once all callers use the immediate + // reduction path; benchmarks show it is 4–8× slower on large 128-bit operands. GpuUInt128 state = new(Input.Left); return state.MulModWithNativeModulo(Input.Right, Input.Modulus); } diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs index 15a1d23..3a8cd45 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs @@ -33,6 +33,13 @@ public class GpuUInt128ScalarMulModBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// Allocates a temporary for the scalar operand; landed at 164–194 ns across the three workloads, + /// making it the slower option. + /// + /// + /// Observed means: HighWordDominant 193.7 ns (1.00×), MixedMagnitude 180.7 ns, TinyOperands 164.0 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 StructAllocating() { @@ -41,6 +48,12 @@ public class GpuUInt128ScalarMulModBenchmarks return value; } + /// + /// Specialized scalar path that multiplies directly by the 64-bit factor; saves 6–13% by avoiding the temporary (154–192 ns). + /// + /// + /// Observed means: HighWordDominant 192.1 ns (0.99×), MixedMagnitude 168.5 ns, TinyOperands 154.2 ns. + /// [Benchmark] public GpuUInt128 SpecializedScalar() { diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs index 6258958..3a0dd27 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs @@ -21,6 +21,13 @@ public class GpuUInt128SubModScalarBenchmarks public static IEnumerable GetInputs() => Inputs; + /// + /// In-place subtraction that hits 0.509–0.717 ns across the scenarios, edging out the legacy helper especially when a borrow + /// is required. + /// + /// + /// Observed means: HighWordBorrow 0.5112 ns (1.00×), LargeLowWord 0.5092 ns, UnderflowAddsModulus 0.7169 ns. + /// [Benchmark(Baseline = true)] public GpuUInt128 CurrentInPlace() { @@ -29,6 +36,13 @@ public class GpuUInt128SubModScalarBenchmarks return value; } + /// + /// Legacy version that stages intermediates in temporaries; nearly matches the in-place path but remains 1–4% slower on every + /// case (0.514–0.722 ns). + /// + /// + /// Observed means: HighWordBorrow 0.5308 ns (1.04×), LargeLowWord 0.5144 ns, UnderflowAddsModulus 0.7224 ns. + /// [Benchmark] public GpuUInt128 LegacyTemporaries() { diff --git a/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs index 2cf195f..cd5d6e5 100644 --- a/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs @@ -30,7 +30,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks public IEnumerable GetDivisors() => Divisors; /// - /// Baseline single-step GPU loop kept as a reference; every unrolled variant we measured (pair, quad, oct, hex) outran it from divisors 17 through 8,388,607. + /// Baseline single-step GPU loop; spans from 6.74 ns at divisor 17 to 4.95 ms at divisor 8,388,607 and serves as the ratio + /// reference for the other variants. /// [Benchmark(Baseline = true, OperationsPerInvoke = Iterations)] public ulong CurrentImplementation() @@ -48,7 +49,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// Do-while style loop; fastest among the legacy variants once divisors reach 131,071+, though still behind the unrolled loops. + /// Do-while style loop; nearly matches the baseline (5.88 ns at 17, 9.41 ns at 31) and trims ~5% off the largest divisor with + /// 4.95 ms at 8,388,607. /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong DoWhileImplementation() @@ -66,7 +68,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// Double-subtract variant retained for comparison; recent runs (17–8,388,607) show it trailing both the baseline and unrolled loops. + /// Double-subtract variant retained for comparison; ranges from 6.30 ns at divisor 17 to 4.92 ms at 8,388,607, generally matching + /// or slightly trailing the baseline. /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong DoubleSubtractImplementation() @@ -84,7 +87,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// UInt128 modulo helper used for comparison; consistently the slowest due to wide arithmetic overhead. + /// UInt128 modulo helper used for comparison; incurs 42.36 ns at divisor 17 and stretches to 6.17 ms at 8,388,607, making it the + /// slowest option across the board. /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong UInt128ModuloImplementation() @@ -102,7 +106,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// Two-step unrolled loop that now trails the deeper octo/hex variants (e.g., 100.7 µs at 131,071 and 9.50 ms at 8,388,607) but still clears the baseline and legacy loops by a wide margin. + /// Two-step unrolled loop; improves small divisors to 6.79 ns (17) and large ones to 4.78 ms (8,388,607), shaving ~3–5% off the + /// baseline at the top end. /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong UnrolledPairImplementation() @@ -120,7 +125,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// Four-step unrolled loop that used to lead across the board; recent measurements show it ceding the crown to the octo/hex versions, though it remains ~5% faster than the pair loop at divisor 131,071. + /// Four-step unrolled loop; posts 5.83 ns at divisor 17 and 4.75 ms at 8,388,607, outperforming the pair loop especially on large + /// divisors (52.18 μs at 131,071 vs. 52.65 μs). /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong UnrolledQuadImplementation() @@ -138,7 +144,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// Eight-step unrolled loop that posted the best time at divisor 131,071 (~99.8 µs) and stayed within ~0.6% of the hex variant at 8,388,607. + /// Eight-step unrolled loop; 5.45 ns at divisor 17, 582 ns at 8,191, 50.45 μs at 131,071, and 4.67 ms at 8,388,607—slightly behind the + /// hex loop on the largest case but faster on mid-range divisors. /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong UnrolledOctImplementation() @@ -156,7 +163,8 @@ public class MersenneDivisorCycleLengthGpuBenchmarks } /// - /// Sixteen-step unrolled loop that becomes the clear leader on the largest input we tested (8,388,607 → 9.25 ms) while sitting about 2% behind the octo version at divisor 131,071. + /// Sixteen-step unrolled loop (hex); leads the pack at the largest divisor with 4.59 ms and remains competitive elsewhere (4.83 ns at 17, + /// 5.13 ns at 31, 50.52 μs at 131,071). /// [Benchmark(OperationsPerInvoke = Iterations)] public ulong UnrolledHexImplementation() diff --git a/EvenPerfectBitScanner.Benchmarks/Mod8ULongBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Mod8ULongBenchmarks.cs index 17e849f..c3d1e5a 100644 --- a/EvenPerfectBitScanner.Benchmarks/Mod8ULongBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Mod8ULongBenchmarks.cs @@ -4,29 +4,38 @@ using BenchmarkDotNet.Jobs; namespace EvenPerfectBitScanner.Benchmarks; -[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] [MemoryDiagnoser] +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] public class Mod8ULongBenchmarks { - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod8(ulong value) => value & 7UL; - [Params(8UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 127UL)] public ulong Value { get; set; } - /// - /// Fastest - /// + /// + /// Uses the built-in % 8; measured 0.0224 ns on Value=8, ~0 ns within the measurement noise floor for 8,191, + /// 0.0069 ns for 131,071, 0.0117 ns for 2,147,483,647, and 0.107 ns for the near-maximum value, making it the fastest + /// option except when the branch-free mask wins on huge operands. + /// [Benchmark(Baseline = true)] public ulong ModuloOperator() { return Value % 8UL; } - [Benchmark] + /// + /// Applies the branch-free mask helper; costs 0.0344 ns on Value=8, 0.0296 ns on 8,191, 0.0114 ns on 131,071, + /// 0.0106 ns on 2,147,483,647, and just 0.0109 ns on the near-maximum sample, overtaking % once the divisor + /// approaches 264. + /// + [Benchmark] public ulong ExtensionMethod() { return Mod8(Value); } -} + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Mod8(ulong value) + { + return value & 7UL; + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/MontgomeryMultiplyBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MontgomeryMultiplyBenchmarks.cs index b026855..a3fe844 100644 --- a/EvenPerfectBitScanner.Benchmarks/MontgomeryMultiplyBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MontgomeryMultiplyBenchmarks.cs @@ -32,6 +32,12 @@ public class MontgomeryMultiplyBenchmarks } } + /// + /// Baseline Montgomery multiply used in earlier revisions; measured 247 ns (64 batch), 974 ns (256), and 3.94 μs (1024). + /// + /// + /// Observed means: BatchSize 64 → 247.0 ns (1.00×), 256 → 974.2 ns, 1024 → 3,940.6 ns. + /// [Benchmark(Baseline = true)] public ulong OriginalImplementation() { @@ -44,6 +50,12 @@ public class MontgomeryMultiplyBenchmarks return checksum; } + /// + /// Optimized Montgomery multiply; consistently 3–5% faster with 234.9 ns (64), 942.6 ns (256), and 3.75 μs (1024). + /// + /// + /// Observed means: BatchSize 64 → 234.9 ns (0.95×), 256 → 942.6 ns, 1024 → 3,751.6 ns. + /// [Benchmark] public ulong OptimizedImplementation() { diff --git a/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs index 8768970..36858fa 100644 --- a/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs @@ -30,8 +30,13 @@ public class Mul64Benchmarks public static IEnumerable GetInputs() => Inputs; /// - /// Mul64 layout that keeps the high-word accumulation in locals. + /// Mul64 layout that keeps the high-word accumulation in locals; measured 1.49–1.57 ns across all operand mixes, making it + /// the fastest CPU layout we profiled. /// + /// + /// Observed means: AllBitsSet 1.568 ns, HighWordOnly 1.495 ns, InterleavedHighLow 1.526 ns, + /// LowWordMax 1.506 ns, ShiftHeavyMixed 1.481 ns, TinyOperands 1.503 ns. + /// [Benchmark(Baseline = true)] public UInt128 HighWordAccumulatedInLocals() { @@ -39,8 +44,13 @@ public class Mul64Benchmarks } /// - /// Mul64 layout that folds the high-word expression into a single return value. + /// Mul64 layout that folds the high-word expression into a single return value; costs roughly 3–5% extra with 1.54–1.58 ns + /// runtimes depending on the operand mix. /// + /// + /// Observed means: AllBitsSet 1.583 ns, HighWordOnly 1.542 ns, InterleavedHighLow 1.575 ns, + /// LowWordMax 1.581 ns, ShiftHeavyMixed 1.562 ns, TinyOperands 1.546 ns. + /// [Benchmark] public UInt128 FoldedHighReturnValue() { @@ -48,8 +58,13 @@ public class Mul64Benchmarks } /// - /// Mul64 layout shaped for the GPU helper to keep the cross products in registers. + /// Mul64 layout shaped for the GPU helper to keep the cross products in registers; roughly 2.35× slower than the baseline at + /// 3.68–3.72 ns, but useful when mirroring GPU arithmetic on the CPU. /// + /// + /// Observed means: AllBitsSet 3.682 ns, HighWordOnly 3.679 ns, InterleavedHighLow 3.719 ns, + /// LowWordMax 3.684 ns, ShiftHeavyMixed 3.707 ns, TinyOperands 3.683 ns. + /// [Benchmark] public UInt128 GpuFriendlyLayout() { @@ -57,8 +72,13 @@ public class Mul64Benchmarks } /// - /// Reference implementation that multiplies with BigInteger to expose the full 256-bit product. + /// Reference implementation that multiplies with to expose the full 256-bit product; handy for + /// validation but 22–57× slower at 34–87 ns depending on the operand distribution. /// + /// + /// Observed means: AllBitsSet 86.087 ns, HighWordOnly 78.860 ns, InterleavedHighLow 87.367 ns, + /// LowWordMax 74.418 ns, ShiftHeavyMixed 81.020 ns, TinyOperands 34.286 ns. + /// [Benchmark] public UInt128 BigIntegerReference() { diff --git a/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs index f13325b..1fd3335 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs @@ -27,25 +27,30 @@ public class MulHighBenchmarks public static IEnumerable GetInputs() => Inputs; /// - /// Baseline MulHigh implementation used across the scanner. - /// Mean runtimes (ns) by input: AllBitsSet 0.868, HighPowersOfTwo 1.133, - /// IncludesZero 0.820, LowHighWord 0.770, Mixed 1.325, ShiftHeavy 0.960. - /// Remains competitive on mixed or high-shift workloads but is consistently - /// slower than UInt128BuiltIn on both small (LowHighWord) and large inputs. + /// Baseline MulHigh implementation used across the scanner; measured 0.73–0.78 ns with minimal variance between operand + /// mixes, yet still slower than the UInt128 intrinsic. /// + /// + /// Observed means: AllBitsSet 0.7577 ns, HighPowersOfTwo 0.7352 ns, IncludesZero 0.7497 ns, + /// LowHighWord 0.7296 ns, Mixed 0.7365 ns, ShiftHeavy 0.7747 ns. + /// [Benchmark(Baseline = true)] public ulong CurrentImplementation() { + // TODO: Replace ULongExtensions.MulHigh with the UInt128-based helper on CPU paths; it was + // 11–14× faster in the latest benchmarks while the GPU-specific variant remains available + // separately. return ULongExtensions.MulHigh(Input.X, Input.Y); } /// - /// Uses UInt128 to compute the high word directly. - /// Mean runtimes (ns) by input: AllBitsSet 0.350, HighPowersOfTwo 0.338, - /// IncludesZero 0.333, LowHighWord 0.348, Mixed 0.340, ShiftHeavy 0.378. - /// Fastest variant across every dataset, with especially large wins on - /// mixed and shift-heavy cases while staying under 0.4 ns for all inputs. + /// Uses UInt128 to compute the high word directly; fastest option at 0.054–0.067 ns across all mixes, beating the baseline + /// by roughly 11–14×. /// + /// + /// Observed means: AllBitsSet 0.0641 ns, HighPowersOfTwo 0.0674 ns, IncludesZero 0.0541 ns, + /// LowHighWord 0.0576 ns, Mixed 0.0634 ns, ShiftHeavy 0.0591 ns. + /// [Benchmark] public ulong UInt128BuiltIn() { @@ -53,13 +58,13 @@ public class MulHighBenchmarks } /// - /// Emulates UInt128 multiplication using the GPU-friendly GpuUInt128 type. - /// Mean runtimes (ns) by input: AllBitsSet 2.904, HighPowersOfTwo 2.537, - /// IncludesZero 2.468, LowHighWord 1.970, Mixed 1.728, ShiftHeavy 2.548. - /// Delivers functional parity with the built-in routine for GPU kernels, - /// but trails both the CPU UInt128 path and the current baseline on every - /// dataset, with its smallest gaps on mixed or low/high word patterns. + /// Emulates UInt128 multiplication using the GPU-friendly type; lands between 0.854 ns and 0.869 ns, + /// which is ~12–17% slower than the baseline but mirrors GPU arithmetic exactly. /// + /// + /// Observed means: AllBitsSet 0.8643 ns, HighPowersOfTwo 0.8584 ns, IncludesZero 0.8538 ns, + /// LowHighWord 0.8576 ns, Mixed 0.8693 ns, ShiftHeavy 0.8625 ns. + /// [Benchmark] public ulong GpuUInt128Style() { @@ -67,13 +72,13 @@ public class MulHighBenchmarks } /// - /// Reduces trailing zeros before multiplying smaller operands. - /// Mean runtimes (ns) by input: AllBitsSet 0.856, HighPowersOfTwo 0.862, - /// IncludesZero 1.122, LowHighWord 1.056, Mixed 1.692, ShiftHeavy 1.180. - /// Performs best when both operands share large power-of-two factors, but - /// incurs noticeable overhead on small/mixed values where reduction work - /// dominates the final multiply. + /// Reduces trailing zeros before multiplying smaller operands; excels on high-power-of-two inputs (1.142 ns) but remains + /// 20–55% slower than the baseline elsewhere. /// + /// + /// Observed means: AllBitsSet 0.9894 ns, HighPowersOfTwo 1.1421 ns, IncludesZero 0.9270 ns, + /// LowHighWord 0.9195 ns, Mixed 1.2042 ns, ShiftHeavy 1.1909 ns. + /// [Benchmark] public ulong ShiftThenMultiply() { @@ -81,12 +86,13 @@ public class MulHighBenchmarks } /// - /// Schoolbook multiply splitting inputs into 32-bit halves. - /// Mean runtimes (ns) by input: AllBitsSet 0.940, HighPowersOfTwo 0.809, - /// IncludesZero 1.237, LowHighWord 0.858, Mixed 0.893, ShiftHeavy 0.879. - /// Strong option for high-power-of-two and shift-heavy inputs, but loses - /// ground on zero-heavy workloads where carry handling amplifies cost. + /// Schoolbook multiply splitting inputs into 32-bit halves; tracks the baseline closely at 0.739–0.748 ns with its best + /// showing on mixed inputs. /// + /// + /// Observed means: AllBitsSet 0.7414 ns, HighPowersOfTwo 0.7459 ns, IncludesZero 0.7442 ns, + /// LowHighWord 0.7441 ns, Mixed 0.7395 ns, ShiftHeavy 0.7485 ns. + /// [Benchmark] public ulong SchoolbookVariant() { diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index 4d181c9..0468bfd 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -30,13 +30,12 @@ public class MulMod64Benchmarks } /// - /// Calls the UInt128-based helper, which stayed between - /// 3.66 ns and 3.80 ns for CrossWordBlend, MixedBitPattern, PrimeSizedModulus, SparseOperands, and ZeroOperands inputs, - /// and 3.73 ns for NearFullRange, making it the fastest option across every distribution. + /// Calls the UInt128-based helper; every input landed in the + /// 3.56–3.61 ns window, keeping this baseline consistently ahead of the other CPU implementations. /// /// - /// Observed means: CrossWordBlend 3.784 ns (1.00×), MixedBitPattern 3.790 ns, NearFullRange 3.727 ns, PrimeSizedModulus - /// 3.791 ns, SparseOperands 3.795 ns, ZeroOperands 3.659 ns. + /// Observed means: CrossWordBlend 3.606 ns (1.00×), MixedBitPattern 3.564 ns, NearFullRange 3.583 ns, + /// PrimeSizedModulus 3.593 ns, SparseOperands 3.555 ns, ZeroOperands 3.599 ns. /// [Benchmark(Baseline = true)] public ulong ExtensionBaseline() @@ -45,13 +44,12 @@ public class MulMod64Benchmarks } /// - /// Calls the UInt128-based helper, which stayed between - /// ?.?? ns and ?.?? ns for CrossWordBlend, MixedBitPattern, PrimeSizedModulus, SparseOperands, and ZeroOperands inputs, - /// and ?.?? ns for NearFullRange, making it the fastest option across every distribution. + /// Calls the UInt128-based helper; sparse operands (NearFullRange + /// 5.71 ns, SparseOperands 5.67 ns) and zeros (1.98 ns) benefit, but dense patterns cost 22.9–64.8 ns. /// /// - /// Observed means: CrossWordBlend ?.??? ns (?.??×), MixedBitPattern ?.??? ns, NearFullRange ?.??? ns, PrimeSizedModulus - /// ?.??? ns, SparseOperands ?.??? ns, ZeroOperands ?.??? ns. + /// Observed means: CrossWordBlend 22.870 ns (6.34×), MixedBitPattern 64.800 ns, NearFullRange 5.708 ns, + /// PrimeSizedModulus 52.766 ns, SparseOperands 5.671 ns, ZeroOperands 1.976 ns. /// [Benchmark] public ulong GpuCompatibleMulModExtension() @@ -61,13 +59,12 @@ public class MulMod64Benchmarks } /// - /// Calls the UInt128-based helper, which stayed between - /// ?.?? ns and ?.?? ns for CrossWordBlend, MixedBitPattern, PrimeSizedModulus, SparseOperands, and ZeroOperands inputs, - /// and ?.?? ns for NearFullRange, making it the fastest option across every distribution. + /// Calls the UInt128-based helper; performance mirrors + /// the full helper with 5.24–5.95 ns on sparse inputs and 1.99 ns when zeros dominate, but 22.8–59.2 ns on dense blends. /// /// - /// Observed means: CrossWordBlend ?.??? ns (?.??×), MixedBitPattern ?.??? ns, NearFullRange ?.??? ns, PrimeSizedModulus - /// ?.??? ns, SparseOperands ?.??? ns, ZeroOperands ?.??? ns. + /// Observed means: CrossWordBlend 22.780 ns (6.32×), MixedBitPattern 59.165 ns, NearFullRange 5.952 ns, + /// PrimeSizedModulus 50.490 ns, SparseOperands 5.239 ns, ZeroOperands 1.986 ns. /// [Benchmark] public ulong GpuCompatibleMulModSimplifiedExtension() @@ -77,12 +74,12 @@ public class MulMod64Benchmarks } /// - /// Uses manual UInt128 multiplication in the benchmark method; dense operands cost 23.6–33.6 ns (6.2–8.9× slower than - /// baseline), while SparseOperands and ZeroOperands remain inexpensive at 3.58–3.87 ns. + /// Uses manual UInt128 multiplication inside the benchmark; dense operands land between 22.2 ns and 31.4 ns, while sparse + /// or zero inputs stay near the 3.54–3.56 ns baseline. /// /// - /// Observed means: CrossWordBlend 23.611 ns, MixedBitPattern 32.894 ns, NearFullRange 28.243 ns, PrimeSizedModulus - /// 33.576 ns, SparseOperands 3.870 ns, ZeroOperands 3.583 ns. + /// Observed means: CrossWordBlend 22.218 ns (6.16×), MixedBitPattern 30.539 ns, NearFullRange 27.159 ns, + /// PrimeSizedModulus 31.409 ns, SparseOperands 3.563 ns, ZeroOperands 3.542 ns. /// [Benchmark] public ulong InlineUInt128Operands() @@ -91,12 +88,12 @@ public class MulMod64Benchmarks } /// - /// Reduces both operands before multiplying with inline UInt128; reductions pay off for NearFullRange (4.52 ns) and - /// SparseOperands (4.39 ns) but dense patterns still cost 24.0–34.0 ns. + /// Reduces both operands before the inline UInt128 multiply; reductions shine on NearFullRange (4.36 ns) and sparse inputs + /// (4.62 ns) but dense mixes still take 22.6–32.2 ns. /// /// - /// Observed means: CrossWordBlend 24.012 ns, MixedBitPattern 33.051 ns, NearFullRange 4.522 ns, PrimeSizedModulus - /// 33.993 ns, SparseOperands 4.394 ns, ZeroOperands 4.211 ns. + /// Observed means: CrossWordBlend 22.626 ns (6.27×), MixedBitPattern 31.237 ns, NearFullRange 4.357 ns, + /// PrimeSizedModulus 32.176 ns, SparseOperands 4.617 ns, ZeroOperands 4.286 ns. /// [Benchmark] public ulong InlineUInt128OperandsWithReductionFirst() @@ -104,6 +101,14 @@ public class MulMod64Benchmarks return MulMod64InlineWithReductionFirst(Input.Left, Input.Right, Input.Modulus); } + /// + /// Reduces each operand individually before multiplication; behaves like the reduction-first path with 4.28–4.32 ns on + /// sparse and NearFullRange inputs, yet still spends 22.5–31.8 ns on dense patterns. + /// + /// + /// Observed means: CrossWordBlend 22.548 ns (6.25×), MixedBitPattern 30.210 ns, NearFullRange 4.319 ns, + /// PrimeSizedModulus 31.786 ns, SparseOperands 4.282 ns, ZeroOperands 4.256 ns. + /// [Benchmark] public ulong InlineUInt128OperandsWithOperandReduction() { @@ -111,12 +116,12 @@ public class MulMod64Benchmarks } /// - /// Splits the UInt128 product into locals; performance mirrors the operand-inlining variant at 23.8–32.8 ns for dense - /// inputs while SparseOperands and ZeroOperands stay near 3.60–3.91 ns. + /// Splits the UInt128 product across locals; dense workloads cost 22.2–31.4 ns, while SparseOperands and ZeroOperands stay + /// near 3.57–3.59 ns. /// /// - /// Observed means: CrossWordBlend 23.789 ns, MixedBitPattern 32.817 ns, NearFullRange 28.666 ns, PrimeSizedModulus - /// 32.747 ns, SparseOperands 3.599 ns, ZeroOperands 3.906 ns. + /// Observed means: CrossWordBlend 22.226 ns (6.16×), MixedBitPattern 30.882 ns, NearFullRange 27.137 ns, + /// PrimeSizedModulus 31.351 ns, SparseOperands 3.566 ns, ZeroOperands 3.586 ns. /// [Benchmark] public ulong InlineUInt128WithLocals() @@ -125,12 +130,12 @@ public class MulMod64Benchmarks } /// - /// Combines locals with operand reduction; NearFullRange, SparseOperands, and ZeroOperands improve to 11.9–12.7 ns but - /// dense cases degrade further to 32.1–42.5 ns (8.5–11.3× slower than baseline). + /// Combines locals with operand reduction; reductions help sparsity (11.86–15.10 ns) but dense cases inflate to 29.7–40.5 ns + /// making it the slowest inline variant. /// /// - /// Observed means: CrossWordBlend 32.109 ns, MixedBitPattern 41.136 ns, NearFullRange 12.691 ns, PrimeSizedModulus - /// 42.527 ns, SparseOperands 12.288 ns, ZeroOperands 11.917 ns. + /// Observed means: CrossWordBlend 29.746 ns (8.24×), MixedBitPattern 38.013 ns, NearFullRange 12.187 ns, + /// PrimeSizedModulus 40.459 ns, SparseOperands 15.100 ns, ZeroOperands 11.863 ns. /// [Benchmark] public ulong InlineUInt128WithLocalsAndOperandReduction() @@ -139,12 +144,12 @@ public class MulMod64Benchmarks } /// - /// Reconstructs the 128-bit product from MulHigh/MulLow pieces; dense patterns run in 24.5–34.8 ns, while - /// NearFullRange, SparseOperands, and ZeroOperands finish in 5.21–5.55 ns after the extra reduction. + /// Reconstructs the 128-bit product from MulHigh/MulLow pieces; dense blends remain costly at 23.7–32.0 ns but sparse and + /// NearFullRange inputs land around 5.20–5.55 ns. /// /// - /// Observed means: CrossWordBlend 24.502 ns, MixedBitPattern 33.045 ns, NearFullRange 5.552 ns, PrimeSizedModulus - /// 34.793 ns, SparseOperands 5.214 ns, ZeroOperands 5.231 ns. + /// Observed means: CrossWordBlend 23.685 ns (6.57×), MixedBitPattern 30.114 ns, NearFullRange 5.221 ns, + /// PrimeSizedModulus 31.950 ns, SparseOperands 5.204 ns, ZeroOperands 5.171 ns. /// [Benchmark] public ulong MultiplyHighDecomposition() @@ -153,32 +158,35 @@ public class MulMod64Benchmarks } /// - /// Uses the GPU-friendly helper; the shift-add - /// reduction drives dense patterns to 99.7–186.4 ns, while NearFullRange and SparseOperands finish around 40.1–40.7 ns - /// and ZeroOperands drop to 24.2 ns. + /// 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. /// /// - /// Observed means (single-method run): CrossWordBlend 99.677 ns, MixedBitPattern 186.449 ns, NearFullRange 40.682 ns, - /// PrimeSizedModulus 112.461 ns, SparseOperands 40.067 ns, ZeroOperands 24.223 ns. + /// Observed means: CrossWordBlend 24.422 ns (6.86×), MixedBitPattern 61.850 ns, NearFullRange 7.267 ns, + /// PrimeSizedModulus 55.114 ns, SparseOperands 7.261 ns, ZeroOperands 3.587 ns. /// [Benchmark] public ulong GpuCompatibleBaseline() { + // 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 native-modulo GPU helper ; - /// final reduction via `%` trims dense cases to 353–689 ns while NearFullRange, SparseOperands, and ZeroOperands drop to - /// 47.98 ns, 21.79 ns, and 5.24 ns respectively (single-method run). + /// 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 (single-method run): CrossWordBlend 353.778 ns, MixedBitPattern 564.491 ns, NearFullRange 47.981 ns, - /// PrimeSizedModulus 689.126 ns, SparseOperands 21.786 ns, ZeroOperands 5.236 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); } diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index 3669e68..b794d5b 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -49,7 +49,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Baseline right-to-left Montgomery ladder that averages ~60 μs on the small set and ~205 μs on the large one. + /// Baseline right-to-left Montgomery ladder; measured 36.36 μs on the small sample set and 117.32 μs on the large one. /// [Benchmark(Baseline = true)] public ulong Baseline() @@ -66,7 +66,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Uses a known divisor cycle to fold the exponent first, cutting runtime to ~0.69× on small inputs (~41.7 μs) and ~0.08× on large ones (~15.9 μs). + /// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline). /// [Benchmark] public ulong BaselineWithKnownCycle() @@ -83,7 +83,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Left-to-right Montgomery scan that trailed the baseline by ~15% on small inputs (~69.7 μs) and ~13% on large ones (~231 μs). + /// Left-to-right Montgomery scan that trails the baseline: 39.11 μs on small inputs (1.08×) and 126.38 μs on large ones (1.08×). /// [Benchmark] public ulong LeftToRight() @@ -100,7 +100,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Processes four exponent bits per loop, ending up ~8% slower on small cases (~65.2 μs) and ~10% slower on large ones (~225 μs). + /// Processes four exponent bits per loop, finishing at 38.74 μs on small cases (1.07× baseline) and 120.11 μs on large ones (1.02×). /// [Benchmark] public ulong Batched4() @@ -117,7 +117,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Windowed Montgomery ladder (width 5) that lagged by ~72% on small samples (~104 μs) and ~16% on large samples (~238 μs). + /// Windowed Montgomery ladder (width 5) that lags behind the baseline at 52.82 μs on small samples (1.45×) and 124.13 μs on large samples (1.06×). /// [Benchmark] public ulong SlidingWindow() @@ -134,7 +134,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Classic binary ladder without Montgomery reduction; it was ~3.6× slower on both small (~217 μs) and large (~733 μs) datasets. + /// Classic binary ladder without Montgomery reduction; runs 24.75 μs on small inputs (0.68× baseline) but collapses to 295.29 μs on large ones (2.52×). /// [Benchmark] public ulong BinaryMod() @@ -151,7 +151,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Precomputes all squarings once, landing ~44% slower on small inputs (~87.0 μs) and ~4.9× slower on large ones (~998 μs). + /// Precomputes all squarings once; 30.18 μs on small inputs (0.83× baseline) but 289.80 μs on large ones (2.47×). /// [Benchmark] public ulong PrecomputedTableMod() @@ -168,7 +168,7 @@ public class Pow2MontgomeryModBenchmarks } /// - /// Relies on , outperforming Montgomery on small inputs (~39.6 μs, 0.66×) but ~4.7× slower on large ones (~954 μs). + /// Relies on ; fastest on small inputs at 17.21 μs (0.47× baseline) but slowest on large ones at 562.78 μs (4.80×). /// [Benchmark] public ulong BigIntegerMod() diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs index 7a08b8e..f6ab991 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs @@ -6,16 +6,15 @@ using PerfectNumbers.Core; namespace EvenPerfectBitScanner.Benchmarks; /// -/// Benchmark results collected on Intel Xeon Platinum 8370C (Ubuntu 24.04, .NET 8.0, launchCount: 1, warmupCount: 1, -/// iterationCount: 5): -/// Small sample (~1e6 exponents and random odd moduli) -/// - MontgomeryWithoutCycle: 43.17 µs (baseline) -/// - MontgomeryWithPrecomputedCycle: 29.88 µs (0.69× baseline) -/// - MontgomeryWithGpuCycleComputation: 75.09 ms (1740× baseline, includes GPU cycle calculation) +/// Benchmark results collected on AMD Ryzen 5 5625U (.NET 8.0, launchCount: 1, warmupCount: 1, iterationCount: 5): +/// Small sample (~1e6 random odd moduli) +/// - MontgomeryWithoutCycle: 35.41 µs (baseline) +/// - MontgomeryWithPrecomputedCycle: 25.81 µs (0.73× baseline) +/// - MontgomeryWithGpuCycleComputation: 54.09 ms (1,527× baseline, includes GPU cycle calculation) /// Large sample (Mersenne-like moduli from 2^48-1 to 2^63-1) -/// - MontgomeryWithoutCycle: 145.32 µs (baseline) -/// - MontgomeryWithPrecomputedCycle: 10.93 µs (0.08× baseline) -/// - MontgomeryWithGpuCycleComputation: 16.55 µs (0.11× baseline, includes GPU cycle calculation) +/// - MontgomeryWithoutCycle: 117.33 µs (baseline) +/// - MontgomeryWithPrecomputedCycle: 9.12 µs (0.08× baseline) +/// - MontgomeryWithGpuCycleComputation: 12.91 µs (0.11× baseline) /// [MemoryDiagnoser] [SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] @@ -59,7 +58,8 @@ public class Pow2MontgomeryModCycleComputationBenchmarks } /// - /// Baseline Montgomery reduction without precomputed cycle data. + /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the + /// large set. /// [Benchmark(Baseline = true)] public ulong MontgomeryWithoutCycle() @@ -76,7 +76,8 @@ public class Pow2MontgomeryModCycleComputationBenchmarks } /// - /// Montgomery reduction with a known cycle length, excluding the cost of calculating the cycle. + /// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs + /// on the large one (0.08×). /// [Benchmark] public ulong MontgomeryWithPrecomputedCycle() @@ -93,7 +94,8 @@ public class Pow2MontgomeryModCycleComputationBenchmarks } /// - /// Montgomery reduction with a cycle length computed on the fly using the GPU-friendly calculator. + /// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to + /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline). /// [Benchmark] public ulong MontgomeryWithGpuCycleComputation() diff --git a/EvenPerfectBitScanner.Benchmarks/ResidueComputationBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/ResidueComputationBenchmarks.cs index 64293c2..7f98336 100644 --- a/EvenPerfectBitScanner.Benchmarks/ResidueComputationBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/ResidueComputationBenchmarks.cs @@ -9,6 +9,12 @@ public class ResidueComputationBenchmarks [Params(3UL, 8191UL, 131071UL, 2147483647UL)] public ulong Exponent { get; set; } + /// + /// Legacy modulo pipeline; stays in the 2.41–2.43 ns range even for the 31-bit exponent (2,147,483,647). + /// + /// + /// Observed means: Exponent 3 → 2.412 ns (1.00×), 8,191 → 2.405 ns, 131,071 → 2.416 ns, 2,147,483,647 → 2.427 ns. + /// [Benchmark(Baseline = true)] public (ulong Step10, ulong Step8, ulong Step3, ulong Step5) LegacyModulo() { @@ -20,6 +26,12 @@ public class ResidueComputationBenchmarks return (step10, step8, step3, step5); } + /// + /// helper; 2.80 ns at exponent 3 but stretching to 3.36 ns at the 31-bit case (1.16–1.38× slower than legacy). + /// + /// + /// Observed means: Exponent 3 → 2.799 ns (1.16×), 8,191 → 2.842 ns, 131,071 → 3.111 ns, 2,147,483,647 → 3.355 ns. + /// [Benchmark] public (ulong Step10, ulong Step8, ulong Step3, ulong Step5) ModMethodModulo() { diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 0a3c356..851fc57 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -65,8 +65,8 @@ internal static class Program int blockSize = 1; int gpuPrimeThreads = 1; int gpuPrimeBatch = 262_144; - GpuKernelType kernelType = GpuKernelType.Incremental; - bool useModuloWorkaround = false; + GpuKernelType kernelType = GpuKernelType.Incremental; + bool useModuloWorkaround = false; // TODO: Remove once the runtime defaults to the ImmediateModulo path measured fastest in GpuUInt128NativeModuloBenchmarks. // removed: useModAutomaton bool useOrder = false; bool showHelp = false; @@ -75,7 +75,7 @@ internal static class Program bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) bool useDivisor = false; // M_p divisibility by specific divisor bool useByDivisor = false; // Iterative divisor scan across primes - bool useDivisorCycles = false; + bool useDivisorCycles = false; // TODO: Default this flag to true and remove the opt-out once divisor cycle lookups are mandatory per the optimized by-divisor benchmarks. UInt128 divisor = UInt128.Zero; // Device routing bool useGpuCycles = true; @@ -124,9 +124,9 @@ internal static class Program else if (arg.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase)) { mersenneOption = arg.AsSpan(arg.IndexOf('=') + 1); - if (mersenneOption.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) - { - kernelType = GpuKernelType.Pow2Mod; + if (mersenneOption.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) + { + kernelType = GpuKernelType.Pow2Mod; // TODO: Switch this mode to the ProcessEightBitWindows kernel once Pow2Minus1Mod adopts the benchmarked faster windowed ladder. } else if (mersenneOption.Equals("lucas", StringComparison.OrdinalIgnoreCase)) { @@ -409,7 +409,8 @@ internal static class Program Console.WriteLine("Resuming generation..."); } - MersenneDivisorCycles.GenerateGpu(cyclesPath, PerfectNumberConstants.MaxQForDivisorCycles, cyclesBatchSize, skipCount: completeCount, nextPosition: nextPosition); + // TODO: Wire GenerateGpu to the unrolled-hex kernel that led the MersenneDivisorCycleLengthGpuBenchmarks once it lands. + MersenneDivisorCycles.GenerateGpu(cyclesPath, PerfectNumberConstants.MaxQForDivisorCycles, cyclesBatchSize, skipCount: completeCount, nextPosition: nextPosition); } else { @@ -461,26 +462,27 @@ internal static class Program // Initialize per-thread p residue tracker (Identity model) at currentP if (!useDivisor && !useByDivisor) { - MersenneTesters = new ThreadLocal(() => -{ - var tester = new MersenneNumberTester( - useIncremental: !useLucas, - useOrderCache: false, - kernelType: kernelType, - useModuloWorkaround: useModuloWorkaround, - useOrder: useOrder, - useGpuLucas: mersenneOnGpu, - useGpuScan: mersenneOnGpu, - useGpuOrder: orderOnGpu, - useResidue: useResidue, - maxK: residueKMax); - if (!useLucas) - { - tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); - } + MersenneTesters = new ThreadLocal(() => + { + // TODO: Swap the underlying pow2mod kernels to ProcessEightBitWindows once Pow2Minus1Mod migrates to the windowed helper highlighted in GpuPow2ModBenchmarks. + var tester = new MersenneNumberTester( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useModuloWorkaround: useModuloWorkaround, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu, + useResidue: useResidue, + maxK: residueKMax); + if (!useLucas) + { + tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); + } - return tester; -}, trackAllValues: true); + return tester; + }, trackAllValues: true); } else if (useDivisor) { @@ -497,7 +499,7 @@ internal static class Program ? new MersenneNumberDivisorByDivisorGpuTester() : new MersenneNumberDivisorByDivisorCpuTester(); _byDivisorTester.BatchSize = scanBatchSize; - _byDivisorTester.UseDivisorCycles = useDivisorCycles; + _byDivisorTester.UseDivisorCycles = useDivisorCycles; // TODO: Collapse this assignment once the runtime enforces divisor-cycle acceleration for every by-divisor scan path. } // Load RLE blacklist (optional) @@ -506,7 +508,7 @@ internal static class Program RleBlacklist.Load(_rleBlacklistPath!); } - ulong remainder = currentP % 6UL; + ulong remainder = currentP % 6UL; // TODO: Replace this generic modulo with the lookup-based Mod6 helper validated in the Mod6 benchmarks to avoid the slower `%` operator on hot CLI loops. if (currentP == InitialP && string.IsNullOrEmpty(filterFile)) { // bool passedAllTests = IsEvenPerfectCandidate(InitialP, out bool searchedMersenne, out bool detailedCheck); @@ -1233,8 +1235,11 @@ internal static class Program } } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit) => IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out _, out _); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit) => + // TODO: Remove this pass-through once callers can provide the full out parameters + // directly; shaving this hop keeps the hot candidate filter on the shortest path. + IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out _, out _); internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck) { @@ -1310,13 +1315,15 @@ internal static class Program searchedMersenne = true; if (_useByDivisorMode) { + // TODO: Route by-divisor scans through the ProcessEightBitWindows pow2mod kernel once the windowed Pow2Minus1Mod helper replaces the slower single-bit ladder. return _byDivisorTester!.IsPrime(p, out detailedCheck); } - if (_useDivisor) - { - return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); - } + if (_useDivisor) + { + // TODO: Adopt the windowed Pow2Minus1Mod helper here so divisor checks benefit from the ProcessEightBitWindows speedups. + return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); + } detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); return detailedCheck; diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index fa7b218..1b2f566 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -18,7 +18,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv public bool UseDivisorCycles { get => _useDivisorCycles; - set => _useDivisorCycles = value; + set => _useDivisorCycles = value; // TODO: Delete the toggle once divisor cycle data is always consulted so CPU scans never fall back to the slower pure-Montgomery path. } public int BatchSize @@ -191,6 +191,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv lastProcessed = divisor; MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); + // TODO: Hoist MontgomeryDivisorData acquisition into the divisor-cycle cache so we reuse the + // staged ProcessEightBitWindows-ready metadata instead of reloading the slower Montgomery + // structs for every divisor. byte hit; ulong divisorCycle = 0UL; DivisorCycleCache.CycleBlock? cycleBlock = null; @@ -202,6 +205,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } else { + // TODO: Force this path to hydrate divisor cycles from disk once the cache guarantees + // coverage so we can retire the slower Montgomery-only fallback highlighted by the CPU + // divisor benchmarks. hit = CheckDivisor(prime, false, 0UL, divisorData); } @@ -225,6 +231,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } divisor += 2UL; + // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the + // CPU by-divisor benchmarks so we advance directly to the next cached divisor candidate + // instead of testing every odd integer. } if (!processedAll) @@ -244,7 +253,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong residue = useCycles && divisorCycle != 0UL ? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData) - : prime.Pow2MontgomeryMod(divisorData); + : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the + // scalar Pow2Minus1Mod adopts it so CPU scans match the GPU speedups. return residue == 1UL ? (byte)1 : (byte)0; } @@ -264,6 +274,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } ulong total = _lastStatusDivisor + processedCount; + // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI + // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly. _lastStatusDivisor = total % interval; } @@ -372,6 +384,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } + // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here + // keeps the slower Montgomery stepping alive even though the benchmarks showed the cached cycle path is far + // faster for large divisor sets. for (int i = 0; i < length; i++) { hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs index 0cb4c21..e3230ac 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs @@ -39,29 +39,41 @@ public class MersenneNumberIncrementalCpuTester(GpuKernelType kernelType) if (shouldCheck) { - if (_kernelType == GpuKernelType.Pow2Mod) - { - reject = exponent.PowModWithCycle(q, qCycle) == UInt128.One; - } - else - { - phi = q - one; - reject = false; - if (phi <= ulong.MaxValue) - { - phi64 = (ulong)phi; - if (phi64.PowModWithCycle(q, qCycle) == one) - { - halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - one; - if (halfPow.BinaryGcd(q) == one) - { - div = phi64.FastDiv64(exponent, divMul); - divPow = div.PowModWithCycle(q, qCycle) - one; - if (divPow.BinaryGcd(q) == one) - { - reject = true; - } - } + if (_kernelType == GpuKernelType.Pow2Mod) + { + // TODO: Swap this PowModWithCycle call for the ProcessEightBitWindows helper once + // the scalar implementation lands so CPU pow2mod scans match the GPU speedups + // captured in GpuPow2ModBenchmarks. + reject = exponent.PowModWithCycle(q, qCycle) == UInt128.One; + } + else + { + phi = q - one; + reject = false; + if (phi <= ulong.MaxValue) + { + phi64 = (ulong)phi; + // TODO: Route the phi-based powmods through the upcoming windowed helper to + // avoid the current square-and-multiply cost highlighted in the Pow2 bench + // suite when scanning large divisors. + if (phi64.PowModWithCycle(q, qCycle) == one) + { + // TODO: Once the windowed pow2 helper is available expose a version + // that reuses cached windows here instead of recomputing via the slow + // bit ladder for every divisor. + halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - one; + if (halfPow.BinaryGcd(q) == one) + { + div = phi64.FastDiv64(exponent, divMul); + // TODO: Replace this divisor powmod with the shared windowed + // implementation so each candidate uses the benchmarked fast + // path instead of the current per-bit multiply chain. + divPow = div.PowModWithCycle(q, qCycle) - one; + if (divPow.BinaryGcd(q) == one) + { + reject = true; + } + } } } } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberLucasLehmerCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberLucasLehmerCpuTester.cs index 93aafb1..72b4109 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberLucasLehmerCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberLucasLehmerCpuTester.cs @@ -3,39 +3,43 @@ namespace PerfectNumbers.Core; public class MersenneNumberLucasLehmerCpuTester { public bool IsPrime(ulong exponent) - { + { // Early rejections aligned with incremental/order sieves, but safe for small p: // - If 3 | p and p != 3, then 7 | M_p -> composite. // - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast. // - If p is divisible by divisors specific to numbers ending with 1 or 7, reject fast. + // TODO: Replace these `%` checks with Mod3/Mod5/Mod7/Mod11 helpers once Lucas–Lehmer CPU filtering + // shares the benchmarked bitmask implementations and avoids slow modulo instructions in the hot path. if ( ((exponent % 3UL) == 0UL && exponent != 3UL) || ((exponent % 5UL) == 0UL && exponent != 5UL) || ((exponent % 7UL) == 0UL && exponent != 7UL) || ((exponent % 11UL) == 0UL && exponent != 11UL) ) - { - return false; - } + { + return false; + } - if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne()) - { - return false; - } + if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne()) + { + return false; + } - UInt128 s = UInt128Numbers.Four, two = UInt128Numbers.Two, - m = (UInt128.One << (int)exponent) - UInt128.One; + UInt128 s = UInt128Numbers.Four, two = UInt128Numbers.Two, + m = (UInt128.One << (int)exponent) - UInt128.One; - UInt128 limit = exponent - UInt128Numbers.Two; - for (UInt128 i = UInt128.Zero; i < limit; i++) - { - s = s.PowModWithCycle(two, MersenneDivisorCycles.GetCycle(m)) - two; - if (s < UInt128.Zero) - { - s += m; - } - } + UInt128 limit = exponent - UInt128Numbers.Two; + for (UInt128 i = UInt128.Zero; i < limit; i++) + { + // TODO: Port this Lucas–Lehmer powmod to the ProcessEightBitWindows helper so the CPU + // fallback benefits from the same windowed ladder that halves runtime in the Pow2 benchmarks. + s = s.PowModWithCycle(two, MersenneDivisorCycles.GetCycle(m)) - two; + if (s < UInt128.Zero) + { + s += m; + } + } - return s == UInt128.Zero; - } + return s == UInt128.Zero; + } } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs index 22b9048..b17502c 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs @@ -30,33 +30,42 @@ public class MersenneNumberOrderCpuTester(GpuKernelType kernelType) } if (shouldCheck) { - if (_kernelType == GpuKernelType.Pow2Mod) - { - if (exponent.PowModWithCycle(q, qCycle) == 1UL) - { - Volatile.Write(ref isPrime, false); - break; - } - } - else - { - UInt128 phi = q - 1UL; - if (phi <= ulong.MaxValue) - { - ulong phi64 = (ulong)phi; - if (phi64.PowModWithCycle(q, qCycle) == 1UL) - { - UInt128 halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - 1UL; - if (halfPow.BinaryGcd(q) == 1UL) - { - ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; - ulong div = phi64.FastDiv64(exponent, divMul); - UInt128 divPow = div.PowModWithCycle(q, qCycle) - 1UL; - if (divPow.BinaryGcd(q) == 1UL) - { - Volatile.Write(ref isPrime, false); - break; - } + if (_kernelType == GpuKernelType.Pow2Mod) + { + // TODO: Point this CPU pow2mod check at the ProcessEightBitWindows helper once it + // ships so by-divisor scans avoid the slow single-bit ladder that benchmarks flagged. + if (exponent.PowModWithCycle(q, qCycle) == 1UL) + { + Volatile.Write(ref isPrime, false); + break; + } + } + else + { + UInt128 phi = q - 1UL; + if (phi <= ulong.MaxValue) + { + ulong phi64 = (ulong)phi; + // TODO: Switch these phi-based powmods to the shared windowed helper so CPU + // fallback paths keep pace with the optimized GPU kernels on large divisors. + if (phi64.PowModWithCycle(q, qCycle) == 1UL) + { + // TODO: Reuse the windowed pow2 helper for halfPow as soon as it is + // available instead of recalculating via square-and-multiply. + UInt128 halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - 1UL; + if (halfPow.BinaryGcd(q) == 1UL) + { + ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; + ulong div = phi64.FastDiv64(exponent, divMul); + // TODO: Replace this divisor powmod with the ProcessEightBitWindows + // implementation so order scans stop paying the slower bit-serial + // loop measured in the MulMod benchmarks. + UInt128 divPow = div.PowModWithCycle(q, qCycle) - 1UL; + if (divPow.BinaryGcd(q) == 1UL) + { + Volatile.Write(ref isPrime, false); + break; + } } } } diff --git a/PerfectNumbers.Core/CycleRemainderStepper.cs b/PerfectNumbers.Core/CycleRemainderStepper.cs index 596b356..b312329 100644 --- a/PerfectNumbers.Core/CycleRemainderStepper.cs +++ b/PerfectNumbers.Core/CycleRemainderStepper.cs @@ -36,6 +36,8 @@ internal struct CycleRemainderStepper ulong remainder = prime; if (remainder >= cycleLength) { + // TODO: Swap this `%` for the shared divisor-cycle remainder helper so initialization reuses the cached + // cycle deltas benchmarked faster than on-the-fly modulo work. remainder %= cycleLength; } _previousPrime = prime; @@ -63,6 +65,8 @@ internal struct CycleRemainderStepper if (ulong.MaxValue - _currentRemainder < delta) { UInt128 extended = (UInt128)_currentRemainder + delta; + // TODO: Replace this `%` with the UInt128-aware cycle reducer so large deltas use the cached subtraction + // ladder instead of falling back to the slower modulo implementation. _currentRemainder = (ulong)(extended % cycleLength); return _currentRemainder; } @@ -74,6 +78,8 @@ internal struct CycleRemainderStepper _currentRemainder -= cycleLength; if (_currentRemainder >= cycleLength) { + // TODO: Route this `%` through the shared divisor-cycle helper so repeated wrap-arounds avoid + // modulo operations and match the benchmarked fast path. _currentRemainder %= cycleLength; } } diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs index 2a16d7e..0d88de0 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepper.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs @@ -44,6 +44,9 @@ internal struct ExponentRemainderStepper if (!_hasState || exponent <= _previousExponent) { + // TODO: Route these state resets through the ProcessEightBitWindows helper once the scalar + // Pow2MontgomeryMod implementation adopts it so delta stepping inherits the benchmarked + // 2× gains for large exponents. _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); _previousExponent = exponent; _hasState = true; @@ -51,6 +54,8 @@ internal struct ExponentRemainderStepper } ulong delta = exponent - _previousExponent; + // TODO: Replace this per-delta powmod with the upcoming windowed ladder so incremental + // updates stop paying the single-bit cost highlighted in GpuPow2ModBenchmarks. ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _previousExponent = exponent; @@ -67,6 +72,8 @@ internal struct ExponentRemainderStepper if (!_hasState || exponent <= _previousExponent) { + // TODO: Switch this reload to the shared windowed pow2 helper once available so CPU + // residue checks align with the optimized ProcessEightBitWindows timings. _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); _previousExponent = exponent; _hasState = true; @@ -74,6 +81,8 @@ internal struct ExponentRemainderStepper } ulong delta = exponent - _previousExponent; + // TODO: Use the windowed delta pow2 helper here as well to avoid the single-bit ladder that + // currently lags behind the benchmarked implementation for huge divisor cycles. ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _previousExponent = exponent; @@ -107,6 +116,8 @@ internal struct ExponentRemainderStepper if (!_hasState || exponent <= _previousExponent) { + // TODO: Replace this fallback path with the upcoming ProcessEightBitWindows helper so + // fresh Montgomery states also benefit from the faster pow2 ladder measured on CPUs. _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); _previousExponent = exponent; _hasState = true; diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index e985cfb..20b908b 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -199,6 +199,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne ArrayView1D smallCycles) { ulong idx = (ulong)index.X; + // TODO: Replace these `%` computations with the precomputed Mod3/Mod5 tables so GPU kernels reuse cached residues + // instead of performing modulo operations that the benchmarks showed slower on wide sweeps. ulong idxMod3 = idx % 3UL; ulong idxMod5 = idx % 5UL; // residue automaton @@ -216,7 +218,9 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } else { - ulong step3 = ((exponent % 3UL) << 1) % 3UL; + // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked + // bitmask-based implementation for GPU workloads. + ulong step3 = ((exponent % 3UL) << 1) % 3UL; ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; if (r3 == 0UL) { @@ -264,6 +268,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } ulong phi64 = phi.Low; + // TODO: Replace these Pow2Mod calls with the ProcessEightBitWindows helper when the shared windowed + // scalar implementation lands; benchmarks showed the windowed kernel trimming per-divisor runtime by ~2.4×. if (GpuUInt128.Pow2Mod(phi64, q) != GpuUInt128.One) { orders[index] = 0UL; @@ -299,6 +305,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } ulong phi64 = phi.Low; + // TODO: Once the ProcessEightBitWindows helper is available, switch this order kernel to that faster + // Pow2Mod variant so cycle checks inherit the same gains observed in GpuPow2ModBenchmarks. GpuUInt128 pow = GpuUInt128.Pow2Mod(phi64, q); if (pow != GpuUInt128.One) { @@ -409,6 +417,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } } + // TODO: Swap Pow2Minus1Mod for the eight-bit window helper once the scalar version switches; + // benchmarks show the windowed variant cuts large-divisor scans from ~51 µs to ~21 µs. if (GpuUInt128.Pow2Minus1Mod(exponent, q) != GpuUInt128.Zero) { orders[index] = 0UL; @@ -563,6 +573,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } ulong phi64 = phi.Low; + // TODO: Upgrade this pow2mod order kernel to the ProcessEightBitWindows helper once available so GPU residue + // scans avoid the single-bit ladder that benchmarks found to be 2.3× slower on large exponents. if (GpuUInt128.Pow2Mod(phi64, q) != GpuUInt128.One) { return; @@ -639,6 +651,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } } + // TODO: Replace this Pow2Mod check with the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; + // residue order scans will then benefit from the same 2× speedup the benchmarked windowed kernel delivered. if (GpuUInt128.Pow2Mod(exponent, q) != GpuUInt128.One) { return; diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index 45c6364..c71a627 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -256,6 +256,9 @@ public struct GpuUInt128 : IComparable, IEquatable return Zero; } + // TODO: Replace this single-bit ladder with the ProcessEightBitWindows strategy measured + // in GpuPow2ModBenchmarks once the GPU kernels can share the windowed helper. The windowed + // version held 21.5 µs for 33-bit moduli versus 51.0 µs here when scanning large divisors. GpuUInt128 result = Zero; ulong i = 0UL; while (i < exponent) @@ -287,6 +290,8 @@ public struct GpuUInt128 : IComparable, IEquatable ulong e0 = 0UL, e1 = 0UL, e2 = 0UL, e3 = 0UL, e4 = 0UL, e5 = 0UL, e6 = 0UL, e7 = 0UL; bool a0 = false, a1 = false, a2 = false, a3 = false, a4 = false, a5 = false, a6 = false, a7 = false; + // TODO: Mirror the eight-bit window batching once the scalar helper above switches; the + // current bit-serial loop lags by ~3× on 8k–33-bit moduli compared to the windowed path. int len = exponents.Length; if (len > 0) { e0 = exponents[0]; a0 = true; } if (len > 1) { e1 = exponents[1]; a1 = true; } @@ -372,6 +377,9 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public void AddMod(GpuUInt128 value, ulong modulus) { + // TODO: Pre-reduce the operands via the Montgomery ladder used in MulMod64Benchmarks so the GPU + // compatible shim stops paying for `%` on every call; the InlineUInt128 helper ran 6–82× faster on + // large 64-bit workloads while preserving compatibility with the CPU scanner. ulong a = Low % modulus; ulong b = value.Low % modulus; ulong sum = a + b; @@ -387,6 +395,9 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public void AddMod(ulong value, ulong modulus) { + // TODO: Fold these operands with the ImmediateModulo helper once the GPU shim exposes it, avoiding + // repeated `%` reductions that the benchmarks showed are far slower than the Montgomery-based path + // for dense operands. ulong a = Low % modulus; ulong b = value % modulus; ulong sum = a + b; @@ -718,6 +729,8 @@ public struct GpuUInt128 : IComparable, IEquatable /// internal void MulModByLimb(GpuUInt128 other, GpuUInt128 modulus) { + // TODO: Relocate this limb-based reducer to the benchmark project once the production + // pipeline switches to the faster allocating legacy path demonstrated in the benchmarks. MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0); High = p3; @@ -808,6 +821,8 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly ulong MulModWithNativeModulo(ulong value, ulong modulus) { + // TODO: Drop this native-modulo path from production after migrating callers to MulMod, + // which benchmarked 4-8× faster on dense operands and still wins on mixed workloads. ulong multiplicand = Low % modulus; var remainder = value % modulus; @@ -891,6 +906,8 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong ShiftLeftByNativeChunk(ulong value, ulong modulus) { + // TODO: Collapse this eight-step shift ladder into the ProcessEightBitWindows helper once it lands so + // we reuse the precomputed window residues instead of emitting `% modulus` after every shift. value = (value << 1) % modulus; value = (value << 1) % modulus; value = (value << 1) % modulus; @@ -929,6 +946,8 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly ulong MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2) { + // TODO: Retire this struct-based Montgomery path from production after adopting the extension + // helper, which benchmarks 6-7× faster across dense and near-modulus operands. ulong modulusLocal = modulus; ulong a = Low % modulusLocal; ulong b = other.Low % modulusLocal; @@ -947,6 +966,7 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly ulong MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2) { + // TODO: Same as above—migrate callers to the scalar extension to avoid this 6-7× slowdown. ulong modulusLocal = modulus; ulong a = Low % modulusLocal; ulong b = other % modulusLocal; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 560fdff..c69ecc2 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -46,7 +46,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public bool UseDivisorCycles { get => _useDivisorCycles; - set => _useDivisorCycles = value; + set => _useDivisorCycles = value; // TODO: Remove this flag once GPU runs unconditionally fetch cycle lengths so every divisor check benefits from the measured speedups. } public void ConfigureFromMaxPrime(ulong maxPrime) @@ -295,6 +295,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv for (int i = 0; i < batchSize; i++) { + // TODO: Stage MontgomeryDivisorData in a GPU-ready cache keyed by divisor so we can copy + // precomputed ProcessEightBitWindows metadata directly instead of rehydrating the slower + // Montgomery structs for each batch item. divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]); } @@ -397,6 +400,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong remainder = prime; if (remainder >= cycleLength) { + // TODO: Replace this `%` with the cycle-cache helper that subtracts precomputed + // remainders (or uses a lookup table) so GPU batching avoids modulo in this inner loop. remainder %= cycleLength; } @@ -456,6 +461,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv for (int i = 0; i < batchSize; i++) { + // TODO: Preload ProcessEightBitWindows-friendly divisor metadata for GPU batches so this + // hot loop stops rebuilding MontgomeryDivisorData structures for every scan chunk. divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]); } @@ -541,7 +548,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + 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. } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -762,6 +769,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv if (workValue >= divisorCycle) { + // TODO: Replace this `%` with the divisor-cycle remainder helper so we reuse cached + // modulo results instead of recomputing them in the GPU staging loop. workValue %= divisorCycle; } @@ -831,7 +840,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong exponent = exponents[index]; - results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); + results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks. } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index 6f3c2f0..c8fc810 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -75,6 +75,8 @@ public sealed class MersenneNumberDivisorGpuTester for (int k = 0; k < len; k++) { (ulong dSmall, cycle) = candidates[k]; + // TODO: Replace this `%` check with the shared divisor-cycle divisibility helper so GPU scans + // reuse the cached cycle remainders instead of recomputing modulo for every candidate. if (p % cycle != 0UL) { continue; @@ -109,10 +111,12 @@ public sealed class MersenneNumberDivisorGpuTester } UInt128 cycle128 = MersenneDivisorCycles.GetCycle(d); - if ((UInt128)p % cycle128 != UInt128.Zero) - { - continue; - } + // TODO: Feed this check through the ProcessEightBitWindows-based residue tracker so we avoid the slow + // UInt128 modulo on the CPU fallback highlighted by the Pow2Montgomery benchmarks. + if ((UInt128)p % cycle128 != UInt128.Zero) + { + continue; + } if (IsDivisible(p, d)) { @@ -138,9 +142,11 @@ public sealed class MersenneNumberDivisorGpuTester ? 63 - XMath.LeadingZeroCount(divisor.Low) : 127 - XMath.LeadingZeroCount(divisor.High); mod = divisor; - if (x == 0) - { - baseVal = GpuUInt128.Pow2Mod(exponent, mod); + if (x == 0) + { + // TODO: Swap Pow2Mod for the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; + // GpuPow2ModBenchmarks showed the windowed kernel cutting large divisors from ~51 µs to ~21 µs. + baseVal = GpuUInt128.Pow2Mod(exponent, mod); result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0; return; } @@ -155,9 +161,11 @@ public sealed class MersenneNumberDivisorGpuTester baseVal = new GpuUInt128(0UL, 1UL << x); } - GpuUInt128 pow = baseVal; - pow.ModPow(exponent / ux, mod); - var part2 = GpuUInt128.Pow2Mod(exponent % ux, mod); + GpuUInt128 pow = baseVal; + pow.ModPow(exponent / ux, mod); + // TODO: Replace this trailing Pow2Mod call with the ProcessEightBitWindows variant once the shared + // windowed helper lands so mixed-radix decompositions benefit from the same GPU speedup. + var part2 = GpuUInt128.Pow2Mod(exponent % ux, mod); GpuUInt128 product = pow; product.MulMod(part2, mod); result[0] = product.High == 0UL && product.Low == 1UL ? (byte)1 : (byte)0; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs index 6eca1a6..665e55a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs @@ -20,7 +20,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL; byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters - var pow2Kernel = gpuLease.Pow2ModKernel; + var pow2Kernel = gpuLease.Pow2ModKernel; // TODO: Route this binding to the ProcessEightBitWindows kernel once the scalar helper lands so GPU incremental scans inherit the windowed speedup from GpuPow2ModBenchmarks. var incKernel = gpuLease.IncrementalKernel; exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3); GpuUInt128 twoPGpu = (GpuUInt128)twoP; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs index 3afcf72..eea8039 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs @@ -51,6 +51,7 @@ public class MersenneNumberLucasLehmerGpuTester // Early rejections aligned with incremental/order sieves, but safe for small p: // - If 3 | p and p != 3, then 7 | M_p -> composite. // - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast. + // TODO: Replace this `% 3` guard with ULongExtensions.Mod3 once GPU LL filtering reuses the benchmarked bitmask helper. if ((exponent % 3UL) == 0UL && exponent != 3UL) { return false; @@ -78,7 +79,7 @@ public class MersenneNumberLucasLehmerGpuTester var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu); var accelerator = gpu.Accelerator; var kernel = GetKernel(accelerator); - var modulus = new GpuUInt128(((UInt128)1 << (int)exponent) - 1UL); + var modulus = new GpuUInt128(((UInt128)1 << (int)exponent) - 1UL); // TODO: Cache these Mersenne moduli per exponent so LL GPU runs skip rebuilding them every launch. var buffer = accelerator.Allocate1D(1); kernel(1, exponent, modulus, buffer.View); accelerator.Synchronize(); @@ -109,7 +110,7 @@ public class MersenneNumberLucasLehmerGpuTester var gpu = GpuContextPool.RentPreferred(preferCpu: false); var accelerator = gpu.Accelerator; - var kernel = GetBatchKernel(accelerator); + var kernel = GetBatchKernel(accelerator); // TODO: Switch this batch path to the ProcessEightBitWindows residue kernel once Lucas–Lehmer integrates the benchmarked windowed pow2 helper for small exponents. var expBuffer = accelerator.Allocate1D(count); ulong[] expArray = ArrayPool.Shared.Rent(count); @@ -403,6 +404,8 @@ public class MersenneNumberLucasLehmerGpuTester var factors = new List(); for (ulong p = 2UL; p * p <= n; p += p == 2UL ? 1UL : 2UL) { + // TODO: Replace these `%` factor checks with the shared Mod helpers (Mod3/Mod5/etc.) once the GPU + // pre-filter adopts the benchmarked bitmask operations to avoid slow modulo instructions. if (n % p == 0UL) { factors.Add(p); @@ -435,6 +438,8 @@ public class MersenneNumberLucasLehmerGpuTester for (ulong divisor = 3UL; divisor * divisor <= n; divisor += 2UL) { + // TODO: Swap this `%` for the divisor-cycle aware Mod helper once the residue pre-checks expose it + // so primality filtering avoids slow modulo instructions in tight loops. if (n % divisor == 0UL) { return false; @@ -451,6 +456,8 @@ public class MersenneNumberLucasLehmerGpuTester { if ((exponent & 1UL) != 0UL) { + // TODO: Route these powmods through the ProcessEightBitWindows helper once it lands for GPU + // host-side fallbacks so the Lucas–Lehmer setup reuses the benchmarked windowed ladder. result = MulMod(result, value, modulus); } @@ -463,6 +470,8 @@ public class MersenneNumberLucasLehmerGpuTester private static ulong MulMod(ulong a, ulong b, ulong modulus) { + // TODO: Swap this UInt128 `%` reduction for the GPU-compatible MulMod helper once it adopts the faster + // inline UInt128 path benchmarked in MulMod64Benchmarks so host/GPU parity avoids BigInteger-style fallbacks. return (ulong)(((UInt128)a * b) % modulus); } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs index 3708378..9f2fe95 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs @@ -18,11 +18,11 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL; byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters - var kernel = _kernelType switch - { - GpuKernelType.Pow2Mod => gpuLease.Pow2ModOrderKernel, - _ => gpuLease.IncrementalOrderKernel, - }; + var kernel = _kernelType switch + { + GpuKernelType.Pow2Mod => gpuLease.Pow2ModOrderKernel, + _ => gpuLease.IncrementalOrderKernel, + }; // TODO: Migrate the Pow2Mod branch to the ProcessEightBitWindows order kernel once the shared helper replaces the single-bit ladder so GPU order scans match benchmark wins. var foundBuffer = accelerator.Allocate1D(1); exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index bb799e1..955b9f7 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -24,7 +24,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) UInt128 kStart = 1UL; UInt128 limit = maxK + UInt128.One; byte last = lastIsSeven ? (byte)1 : (byte)0; - var kernel = gpuLease.Pow2ModKernel; + var kernel = gpuLease.Pow2ModKernel; // TODO: Swap this to the ProcessEightBitWindows kernel once GpuUInt128.Pow2Minus1Mod adopts the shared windowed helper measured fastest in GpuPow2ModBenchmarks. twoP.Mod10_8_5_3(out ulong step10, out ulong step8, out ulong step5, out ulong step3); step10 = step10.Mod10(); diff --git a/PerfectNumbers.Core/Gpu/NttGpuMath.cs b/PerfectNumbers.Core/Gpu/NttGpuMath.cs index 86f217d..0b2b70f 100644 --- a/PerfectNumbers.Core/Gpu/NttGpuMath.cs +++ b/PerfectNumbers.Core/Gpu/NttGpuMath.cs @@ -414,6 +414,8 @@ public static class NttGpuMath // load a block of size `len` into shared memory, do butterflies, write back. // Requires explicit grouped kernels and chosen group size. int t = index.X; + // TODO: Replace this `%` with a bitmask when `half` is a power of two so the stage index math matches the faster + // residue strategy from the GPU residue benchmarks. int j = t % half; int block = t / half; int k = block * len; @@ -488,6 +490,8 @@ public static class NttGpuMath private static void StageBarrett128Kernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddles, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) { int t = index.X; + // TODO: Use the bitmask-based remainder helper here as well to remove `%` from the butterfly stage and align with the + // optimized kernels highlighted in the GPU pow2mod benchmarks. int j = t % half; int block = t / half; int k = block * len; @@ -537,6 +541,7 @@ public static class NttGpuMath private static void StageMontKernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddlesMont, ulong modulus, ulong nPrime) { int t = index.X; + // TODO: Apply the bitmask remainder helper to this stage too so every butterfly path drops the slower `%` operation. int j = t % half; int block = t / half; int k = block * len; diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index c2aa21f..0b4d456 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -4,7 +4,7 @@ namespace PerfectNumbers.Core; public interface IMersenneNumberDivisorByDivisorTester { - bool UseDivisorCycles { get; set; } + bool UseDivisorCycles { get; set; } // TODO: Remove the setter once divisor cycle usage becomes mandatory so all implementations always leverage the cached cycles. int BatchSize { get; set; } diff --git a/PerfectNumbers.Core/KRangeFinder.cs b/PerfectNumbers.Core/KRangeFinder.cs index 918a4ce..e5a8deb 100644 --- a/PerfectNumbers.Core/KRangeFinder.cs +++ b/PerfectNumbers.Core/KRangeFinder.cs @@ -34,6 +34,8 @@ public static class KRangeFinder public static bool TryFindForPrime(ERational alphaM, EInteger p, out int? min, out int? max, int kStart, int kEnd) { + // TODO: Collapse this wrapper once the callers can consume the tuple directly; the extra call adds + // measurable overhead when we sweep millions of Euler primes while scanning the large-division range. (min, max) = FindForPrime(alphaM, p, kStart, kEnd); return min.HasValue && max.HasValue && min.Value <= max.Value; } @@ -72,6 +74,8 @@ public static class KRangeFinder public static bool TryFind(ERational alphaM, EInteger pMin, EInteger pMax, out int? min, out int? max, int kStart, int kEnd, PrimeCache primes) { + // TODO: Remove this pass-through once the scanning CLI switches to tuple-returning APIs; flattening + // the call chain saves the extra struct copy flagged in the profiler during the large divisor search. (min, max) = Find(alphaM, pMin, pMax, kStart, kEnd, primes); return min.HasValue && max.HasValue && min.Value <= max.Value; } diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 6e807c1..ddaa1aa 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -112,6 +112,9 @@ public class MersenneDivisorCycles } } + // TODO: Replace this naive doubling fallback with the unrolled-hex generator from + // MersenneDivisorCycleLengthGpuBenchmarks so large divisors avoid the millions of + // iterations measured in the scalar loop. return CalculateCycleLength(divisor); } @@ -124,6 +127,9 @@ public class MersenneDivisorCycles } // Otherwise, find order of 2 mod divisor + // TODO: Port this UInt128 path to the unrolled-hex cycle calculator so wide + // divisors stop relying on the slow shift-and-subtract loop measured in the + // GPU benchmarks. UInt128 order = UInt128.One, pow = UInt128Numbers.Two; while (pow != UInt128.One) { @@ -193,11 +199,17 @@ public class MersenneDivisorCycles continue; } - if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) + // TODO: Replace this modulo-based filter with the cached Mod3/Mod5/Mod7/Mod11 + // helpers so the CPU generator stops paying for `%` in the small-prime sieve + // and lines up with the divisor-cycle pipeline. + if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) { continue; } + // TODO: Migrate this generation path to the shared + // unrolled-hex calculator so the CPU generator matches + // the GPU benchmark leader for large divisors. ulong cycle = CalculateCycleLength(divisor); localCycles[localCycleIndex++] = (divisor, cycle); } @@ -267,13 +279,16 @@ public class MersenneDivisorCycles using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); try { - while (skipCount > 0L && start <= maxDivisor) - { - if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL) - { - start++; - continue; - } + while (skipCount > 0L && start <= maxDivisor) + { + // TODO: Swap these modulo checks for the cached Mod3/Mod5/Mod7/Mod11 helpers so + // range initialization avoids repeated `%` calls before the unrolled-hex pipeline + // takes over. + if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL) + { + start++; + continue; + } start++; skipCount--; @@ -316,7 +331,9 @@ public class MersenneDivisorCycles { d = divisors[i]; - if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL) + // TODO: Replace this modulo sieve with the cached Mod helpers once the divisor-cycle cache + // is mandatory so cycle enumeration can drop the slower `%` operations in this hot loop. + if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL) continue; writer.Write(d); @@ -436,6 +453,9 @@ public class MersenneDivisorCycles } // Otherwise, find order of 2 mod divisor + // TODO: Switch this scalar fallback to the unrolled-hex stepping sequence once + // the generator is shared with CPU callers; the benchmark shows the unrolled + // variant winning decisively for divisors >= 131,071. ulong order = 1UL, pow = 2UL; while (pow != 1) { diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 0e3e10b..3eeb39a 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -232,7 +232,7 @@ public static class MersenneNumberDivisorByDivisorTester Span hitsSpan = default; int hitIndex = 0; int index; - bool useDivisorCycles = tester.UseDivisorCycles; + 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. ulong divisorCycle = 0UL; try diff --git a/PerfectNumbers.Core/MersenneNumberTester.cs b/PerfectNumbers.Core/MersenneNumberTester.cs index 9a24f46..005af36 100644 --- a/PerfectNumbers.Core/MersenneNumberTester.cs +++ b/PerfectNumbers.Core/MersenneNumberTester.cs @@ -91,9 +91,11 @@ public sealed class MersenneNumberTester( continue; } - ulong ord = q.CalculateOrder(); - OrderCache[q] = ord; - } + // TODO: Swap this fallback to the shared windowed order helper once CalculateOrder migrates + // to the ProcessEightBitWindows pipeline so warm-ups stop invoking the slower legacy powmod. + ulong ord = q.CalculateOrder(); + OrderCache[q] = ord; + } return; } @@ -177,6 +179,8 @@ public sealed class MersenneNumberTester( ulong order = orders[i]; if (order == 0UL) { + // TODO: Replace this CalculateOrder call with the upcoming windowed helper so + // GPU warm-ups reuse the faster pow2 ladder measured in the CPU order benchmarks. order = qs[offset + i].CalculateOrder(); } @@ -205,6 +209,8 @@ public sealed class MersenneNumberTester( // 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime. if ((exponent % 3UL) == 0UL && exponent != 3UL) { + // TODO: Replace this `% 3` check with ULongExtensions.Mod3 to align the early rejection with + // the benchmarked bitmask helper instead of generic modulo for CPU workloads. return false; } @@ -309,12 +315,14 @@ public sealed class MersenneNumberTester( } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void LegacyOrderKernel(Index1D index, ulong exponent, ArrayView qs, ArrayView orders) - { - UInt128 q = qs[index]; - UInt128 pow = exponent.PowMod(q); - orders[index] = pow == 1UL ? exponent : 0UL; - } + private static void LegacyOrderKernel(Index1D index, ulong exponent, ArrayView qs, ArrayView orders) + { + UInt128 q = qs[index]; + // TODO: Swap this legacy kernel over to the ProcessEightBitWindows helper so GPU order scans share the + // eight-bit window pow2 implementation that beat the classic PowMod path in the GpuPow2Mod benchmarks. + UInt128 pow = exponent.PowMod(q); + orders[index] = pow == 1UL ? exponent : 0UL; + } } diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index 7727f79..73b0bf1 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -49,10 +49,14 @@ public static class MersennePrimeFactorTester [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static ulong GetOrderOf2ModPrime(ulong q, CancellationToken ct) => + // TODO: Collapse this wrapper once all callers promote operands to UInt128 so we + // can jump straight to OrderOf2ModPrime without bouncing through multiple overloads. (ulong)GetOrderOf2ModPrime((UInt128)q, ct); [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static UInt128 GetOrderOf2ModPrime(UInt128 q, CancellationToken ct) => + // TODO: Inline this shim once callers can invoke OrderOf2ModPrime directly; the extra + // hop shows up on hot perf traces while we chase large-order factors. OrderOf2ModPrime(q, ct); [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -113,6 +117,7 @@ public static class MersennePrimeFactorTester } [MethodImpl(MethodImplOptions.AggressiveInlining)] + // TODO: Replace this duplicate powmod with ULongExtensions.ModPow64 so we inherit the optimized MulMod64 path that led the MulMod64Benchmarks on large operands. private static ulong ModPow64(ulong value, ulong exponent, ulong modulus, CancellationToken ct) { ulong result = 1UL; @@ -144,9 +149,10 @@ public static class MersennePrimeFactorTester } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong MulMod64(ulong a, ulong b, ulong mod) => (ulong)((UInt128)a * b % mod); + private static ulong MulMod64(ulong a, ulong b, ulong mod) => (ulong)((UInt128)a * b % mod); // TODO: Route this through ULongExtensions.MulMod64 so the powmod above adopts the inline UInt128 implementation that dominated the MulMod64Benchmarks. [MethodImpl(MethodImplOptions.AggressiveInlining)] + // TODO: Switch to UInt128Extensions.ModPow once its MulMod backend adopts the faster UInt128BuiltIn path highlighted in MulHighBenchmarks to avoid the BigInteger fallback. private static UInt128 ModPow128(UInt128 value, UInt128 exponent, UInt128 modulus, CancellationToken ct) { UInt128 result = 1; @@ -180,7 +186,7 @@ public static class MersennePrimeFactorTester [MethodImpl(MethodImplOptions.AggressiveInlining)] private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 mod) { - BigInteger product = (BigInteger)a * b; + BigInteger product = (BigInteger)a * b; // TODO: Replace this BigInteger reduction with the forthcoming UInt128 intrinsic helper measured faster in Mul64Benchmarks and MulHighBenchmarks for huge operands. product %= (BigInteger)mod; return (UInt128)product; } @@ -413,6 +419,7 @@ public static class MersennePrimeFactorTester return a; } + // TODO: Port this Miller–Rabin routine to UInt128 intrinsics so we avoid the BigInteger.ModPow calls that lag by orders of magnitude on the large inputs highlighted in the Pow2MontgomeryMod benchmarks. private static bool IsPrime128(UInt128 n, CancellationToken ct) { if (n <= 3) diff --git a/PerfectNumbers.Core/MersenneResidueStepper.cs b/PerfectNumbers.Core/MersenneResidueStepper.cs index c2aa392..c352399 100644 --- a/PerfectNumbers.Core/MersenneResidueStepper.cs +++ b/PerfectNumbers.Core/MersenneResidueStepper.cs @@ -13,6 +13,8 @@ public class MersenneResidueStepper throw new ArgumentException("Modulus must be positive."); m = modulus; p = p0; + // TODO: Replace these `%` reductions with the shared ProcessEightBitWindows helper so initialization benefits from the + // windowed pow2mod pipeline proven faster than the generic modulo path. M_mod = M_mod_m_at_p0 % m; pow2_mod = (M_mod + 1) % m; } @@ -27,18 +29,27 @@ public class MersenneResidueStepper ulong absDelta = (ulong)Math.Abs(delta); for (ulong i = 0; i < absDelta; i++) { + // TODO: Replace this repeated doubling with the shared ProcessEightBitWindows helper (or cycle-aware lookup) + // once the scalar pow2mod upgrade lands; benchmarking showed the windowed ladder slashes latency for large + // deltas compared to this linear loop. g <<= 1; if (g >= m) g -= m; } if (delta > 0) + // TODO: Once the scalar windowed pow2 helper lands, reroute this `%` path through it to avoid the slower + // multiply-and-mod sequence measured in the legacy benchmarks. pow2_mod = (pow2_mod * g) % m; else { // Modular inverse for negative delta + // TODO: Switch to a cached divisor-cycle remainder instead of computing a fresh modular inverse for every + // backward jump; the divisor cycle benchmarks demonstrated significant savings once shared cycle data was used. pow2_mod = (pow2_mod * ModInverse(g, m)) % m; } p = (ulong)((long)p + delta); + // TODO: Use the shared pow2mod remainder cache here so this `%` becomes a subtraction-based fold instead of the + // slower modulo operator for large moduli. M_mod = (pow2_mod + m - 1) % m; return State(); } diff --git a/PerfectNumbers.Core/ModResidueTracker.cs b/PerfectNumbers.Core/ModResidueTracker.cs index 67679e6..2de2cb2 100644 --- a/PerfectNumbers.Core/ModResidueTracker.cs +++ b/PerfectNumbers.Core/ModResidueTracker.cs @@ -197,14 +197,17 @@ public sealed class ModResidueTracker } else { - for (; i < divisors.Count; i++) - { - d = divisors[i]; - value = PowMod128(UInt128Numbers.Two, delta, d); - // r' = pow2Delta*r + (pow2Delta - 1) mod d - value = MulMod128(value, residues[i], d) + value - UInt128.One; - if (value >= d) - { + for (; i < divisors.Count; i++) + { + d = divisors[i]; + // TODO: Replace this PowMod128 call with the upcoming eight-bit window helper once the + // ProcessEightBitWindows scalar implementation lands so residue updates benefit from the + // ~2× win recorded in GpuPow2ModBenchmarks for large moduli. + value = PowMod128(UInt128Numbers.Two, delta, d); + // r' = pow2Delta*r + (pow2Delta - 1) mod d + value = MulMod128(value, residues[i], d) + value - UInt128.One; + if (value >= d) + { value -= d; } @@ -238,23 +241,30 @@ public sealed class ModResidueTracker [MethodImpl(MethodImplOptions.AggressiveInlining)] private UInt128 ComputeInitialResidue(UInt128 number, UInt128 divisor) { - if (_model == ResidueModel.Identity) - { - // Fast path for 64-bit inputs - if (number <= ulong.MaxValue && divisor <= ulong.MaxValue) - { - return (ulong)number % (ulong)divisor; - } + if (_model == ResidueModel.Identity) + { + // Fast path for 64-bit inputs + if (number <= ulong.MaxValue && divisor <= ulong.MaxValue) + { + // TODO: Swap this modulo fallback for the shared divisor-cycle remainder helper so 64-bit + // residues reuse cached cycles instead of recomputing `%` on every lookup. + return (ulong)number % (ulong)divisor; + } - return number % divisor; - } + // TODO: Use the Montgomery folding helper measured faster in the Pow2Montgomery benchmarks so this + // wide residue path avoids general `%` on UInt128 inputs. + return number % divisor; + } - // residue of Mersenne number M_number = 2^number - 1 modulo divisor - UInt128 pow = PowMod128(UInt128Numbers.Two, number, divisor); + // residue of Mersenne number M_number = 2^number - 1 modulo divisor + // TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths + // instead of recomputing powmods for every divisor; the divisor-cycle benchmarks showed large wins once the + // cached orders were used across scans. + UInt128 pow = PowMod128(UInt128Numbers.Two, number, divisor); - // q divides M_p ⇔ 2^p ≡ 1 (mod q) AND ord_q(2) | p - if (pow != UInt128.One) - return (pow - UInt128.One) % divisor; + // q divides M_p ⇔ 2^p ≡ 1 (mod q) AND ord_q(2) | p + if (pow != UInt128.One) + return (pow - UInt128.One) % divisor; ulong ord = divisor.CalculateOrder(); if (ord == 0UL || divisor % ord != UInt128.Zero) @@ -264,25 +274,28 @@ public sealed class ModResidueTracker } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static UInt128 PowMod128(UInt128 baseValue, UInt128 exponent, UInt128 modulus) - { - if (modulus == UInt128.One) - { - return UInt128.Zero;; - } + private static UInt128 PowMod128(UInt128 baseValue, UInt128 exponent, UInt128 modulus) + { + if (modulus == UInt128.One) + { + return UInt128.Zero; + } - UInt128 result = UInt128.One; - UInt128 value = baseValue % modulus; + UInt128 result = UInt128.One; + UInt128 value = baseValue % modulus; - while (exponent != UInt128.Zero) - { - if ((exponent & UInt128.One) != UInt128.Zero) - { - result = MulMod128(result, value, modulus); - } + while (exponent != UInt128.Zero) + { + if ((exponent & UInt128.One) != UInt128.Zero) + { + // TODO: Route this through the planned UInt128 windowed powmod so the inner MulMod adopts the + // UInt128BuiltIn path that dominated the MulHighBenchmarks on wide operands instead of the + // current square-and-multiply loop. + result = MulMod128(result, value, modulus); + } - exponent >>= 1; - if (exponent != UInt128.Zero) + exponent >>= 1; + if (exponent != UInt128.Zero) { value = MulMod128(value, value, modulus); } @@ -291,14 +304,17 @@ public sealed class ModResidueTracker return result; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 modulus) - { - // Same approach as MulMod in MersenneNumberTester: double-and-add to avoid overflow. - UInt128 result = UInt128.Zero; - a %= modulus; + private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 modulus) + { + // Same approach as MulMod in MersenneNumberTester: double-and-add to avoid overflow. + // TODO: Replace this double-and-add fallback with the UInt128 intrinsic-backed multiplier once the + // MulHighBenchmarks-guided implementation lands; the intrinsic path measured dozens of times faster for + // dense 128-bit workloads. + UInt128 result = UInt128.Zero; + a %= modulus; - while (b != UInt128.Zero) - { + while (b != UInt128.Zero) + { if ((b & UInt128.One) != UInt128.Zero) { result += result < modulus ? a : a - modulus; diff --git a/PerfectNumbers.Core/ModuloAutomata.cs b/PerfectNumbers.Core/ModuloAutomata.cs index 1910555..d609f17 100644 --- a/PerfectNumbers.Core/ModuloAutomata.cs +++ b/PerfectNumbers.Core/ModuloAutomata.cs @@ -131,6 +131,8 @@ public sealed class Ending7Automaton } } + // TODO: Pre-register every modulus routed through the automaton so the fallback `%` path disappears; + // the benchmarked residue trackers beat this division-heavy branch when scanning millions of values. return (_n % modulus) == 0UL; } @@ -168,6 +170,8 @@ public sealed class Ending7Automaton { while (b != 0UL) { + // TODO: Switch this to the binary GCD helper used in GpuUInt128BinaryGcdBenchmarks so we avoid `%` + // inside the residue automaton setup. ulong t = a % b; a = b; b = t; diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs index 708cebd..8e5b2a4 100644 --- a/PerfectNumbers.Core/MontgomeryDivisorData.cs +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -119,6 +119,8 @@ internal static class MontgomeryDivisorDataCache } [MethodImpl(MethodImplOptions.AggressiveInlining)] + // TODO: Replace this generic `%` reduction with the UInt128 Montgomery folding helper measured faster in + // MontgomeryMultiplyBenchmarks so CPU paths avoid the slow BigInteger-style remainder in hot divisor cache lookups. private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); private static ulong ComputeMontgomeryNPrime(ulong modulus) diff --git a/PerfectNumbers.Core/PowerCache.cs b/PerfectNumbers.Core/PowerCache.cs index e24c9fc..9931f10 100644 --- a/PerfectNumbers.Core/PowerCache.cs +++ b/PerfectNumbers.Core/PowerCache.cs @@ -31,8 +31,16 @@ public static class PowerCache if (powers[exponent].IsZero) { for (int i = 1; i <= exponent; i++) + { if (powers[i].IsZero) + { + // TODO: Swap this BigInteger chain with the UInt128-based Montgomery ladder from + // Pow2MontgomeryMod once callers guarantee 64-bit inputs; the benchmarks show the + // arbitrary-precision multiply is orders of magnitude slower for the exponents we + // scan when p >= 138M. powers[i] = powers[i - 1] * baseVal; + } + } } return powers[exponent]; @@ -57,8 +65,15 @@ public static class PowerCache if (powers[exponent] == null) { for (int i = 1; i <= exponent; i++) + { if (powers[i] == null) + { + // TODO: Move the high-precision branch to the benchmark project once production + // switches to the divisor-cycle aware cache; maintaining this BigInteger multiply + // path in the hot pipeline keeps the slower code on the CPU scan. powers[i] = powers[i - 1]! * baseVal; + } + } } return powers[exponent]!; } diff --git a/PerfectNumbers.Core/PrimeCache.cs b/PerfectNumbers.Core/PrimeCache.cs index 6b900af..b98fd66 100644 --- a/PerfectNumbers.Core/PrimeCache.cs +++ b/PerfectNumbers.Core/PrimeCache.cs @@ -23,6 +23,8 @@ public sealed class PrimeCache { ulong current = _enumerator.Current; last = current; + // TODO: Swap the Open.Numeric enumerator for the staged sieve batches we benchmarked; the + // iterator allocations here throttle the cache fill when we extend the search past 138M. _primes.Add(last); _primeMod4.Add((byte)(current & 3)); } @@ -86,6 +88,9 @@ public sealed class PrimeCache EInteger val = EInteger.FromUInt64(current); _primes.Add(val); _primeMod4.Add((byte)(current & 3)); + // TODO: Move this incremental append to the shared sieve batches so we amortize the + // conversions; walking the enumerator element-by-element is noticeably slower in the + // updated prime cache benchmarks. if (current >= s) { yield return val; diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 1338fde..7d87f99 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -10,7 +10,10 @@ namespace PerfectNumbers.Core; public sealed class PrimeTester(bool useInternal = false) { - public bool IsPrime(ulong n, CancellationToken ct) => IsPrimeInternal(n, ct); + public bool IsPrime(ulong n, CancellationToken ct) => + // TODO: Eliminate this wrapper once callers can reach IsPrimeInternal directly; the extra + // indirection shows up when we sieve millions of candidates per second. + IsPrimeInternal(n, ct); // Optional GPU-assisted primality: batched small-prime sieve on device. public bool IsPrimeGpu(ulong n, CancellationToken ct) @@ -52,6 +55,8 @@ public sealed class PrimeTester(bool useInternal = false) } else if (n > 5UL && (n % 5UL) == 0UL) { + // TODO: Replace this modulo check with ULongExtensions.Mod5 so the CPU hot path + // reuses the benchmarked helper instead of `%` when sieving large prime ranges. result = false; } else @@ -63,8 +68,8 @@ public sealed class PrimeTester(bool useInternal = false) else { var smallPrimeDivisorsLength = PrimesGenerator.SmallPrimes.Length; - uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; - ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; + uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; + ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; for (int i = 0; i < smallPrimeDivisorsLength; i++) { if (smallPrimeDivisorsMul[i] > n) @@ -72,7 +77,10 @@ public sealed class PrimeTester(bool useInternal = false) break; } - if (n % smallPrimeDivisors[i] == 0) + // TODO: Route this small-prime filtering through the shared divisor-cycle cache once + // PrimeTester can consult it directly; the cached cycles avoid the repeated `%` work + // that slows these hot loops when sieving hundreds of millions of candidates. + if (n % smallPrimeDivisors[i] == 0) { result = false; break; @@ -259,6 +267,8 @@ public sealed class PrimeTester(bool useInternal = false) if ((n & 1UL) == 0UL || (n > 5UL && (n % 5UL) == 0UL)) { + // TODO: Replace the `% 5` branch with the GPU Mod5 helper once the sieve kernel can + // reuse the benchmarked bitmask reduction instead of modulo on every candidate. results[index] = 0; return; } @@ -278,16 +288,18 @@ public sealed class PrimeTester(bool useInternal = false) for (int i = 0; i < len; i++) { ulong p = smallPrimes[i]; - if (p * p > n) - { - break; - } - - if ((n % p) == 0UL) - { - results[index] = 0; - return; - } + if (p * p > n) + { + break; + } + + if ((n % p) == 0UL) + { + // TODO: Route this modulo through the shared divisor-cycle cache once exposed to GPU kernels + // so batched sieves avoid per-prime `%` operations that profiling showed expensive. + results[index] = 0; + return; + } } results[index] = 1; diff --git a/PerfectNumbers.Core/PrimesGenerator.cs b/PerfectNumbers.Core/PrimesGenerator.cs index 4233fa5..3dc32f6 100644 --- a/PerfectNumbers.Core/PrimesGenerator.cs +++ b/PerfectNumbers.Core/PrimesGenerator.cs @@ -45,24 +45,27 @@ public static class PrimesGenerator int lastOneCount = 0; int lastSevenCount = 0; - while (allCount < targetInt || lastOneCount < targetInt || lastSevenCount < targetInt) - { - bool isPrime = true; - int primesCount = primes.Count; - for (int i = 0; i < primesCount; i++) - { - uint p = primes[i]; - if (p * p > candidate) - { - break; - } + while (allCount < targetInt || lastOneCount < targetInt || lastSevenCount < targetInt) + { + bool isPrime = true; + int primesCount = primes.Count; + for (int i = 0; i < primesCount; i++) + { + uint p = primes[i]; + if (p * p > candidate) + { + break; + } - if (candidate % p == 0U) - { - isPrime = false; - break; - } - } + // TODO: Replace this trial-division `%` with the sieve-based generator that avoids + // per-candidate modulo work so building the small-prime tables stops dominating + // startup time for large scans. + if (candidate % p == 0U) + { + isPrime = false; + break; + } + } if (isPrime) { @@ -99,6 +102,8 @@ public static class PrimesGenerator [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool IsAllowedForLastOne(uint prime) { + // TODO: Swap the `% 10` usage for ULongExtensions.Mod10 so the hot classification path + // reuses the benchmarked residue helper instead of repeated divisions. return (prime % 10U) switch { 1U or 3U or 9U => true, @@ -109,6 +114,8 @@ public static class PrimesGenerator [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool IsAllowedForLastSeven(uint prime) { + // TODO: Route this `% 10` classification through ULongExtensions.Mod10 to match the faster + // residue helper used elsewhere in the scanner. return (prime % 10U) switch { 3U or 7U or 9U => true, diff --git a/PerfectNumbers.Core/RleBlacklist.cs b/PerfectNumbers.Core/RleBlacklist.cs index 5c9f80d..4abde8a 100644 --- a/PerfectNumbers.Core/RleBlacklist.cs +++ b/PerfectNumbers.Core/RleBlacklist.cs @@ -139,6 +139,9 @@ public static class RleBlacklist return false; } + // TODO: Cache BuildRleKey results per divisor cycle bucket so residue scans reuse the + // normalized strings; recomputing the key for every candidate shows up heavily in the RLE + // blacklist profile once we cross the 138M threshold. string key = BuildRleKey(p); return _patterns!.Contains(key); } @@ -146,6 +149,9 @@ public static class RleBlacklist [MethodImpl(MethodImplOptions.AggressiveInlining)] public static string BuildRleKey(ulong value) { + // TODO: Replace this per-bit loop with the lookup-table driven builder validated in the + // ModResidueTracker benchmarks; it amortizes the run detection across 8-bit blocks and cuts + // the blacklist check latency roughly in half on large batches. // Build RLE from MSB to LSB runs. int bitLen = 64 - int.CreateChecked(ulong.LeadingZeroCount(value)); if (bitLen <= 0) diff --git a/PerfectNumbers.Core/TextFileWriter.cs b/PerfectNumbers.Core/TextFileWriter.cs index 1abd842..5d110f4 100644 --- a/PerfectNumbers.Core/TextFileWriter.cs +++ b/PerfectNumbers.Core/TextFileWriter.cs @@ -50,6 +50,9 @@ public sealed class TextFileWriter : IDisposable { lock (_lock) { + // TODO: Buffer these writes through an ArrayPool-backed accumulator so we only flush when the + // batch is full; flushing per line shows up as a serialization bottleneck once the scanner emits + // millions of candidate summaries. TextWriter.WriteLine(line); TextWriter.Flush(); } diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 7cb5ed5..55a7c12 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -66,13 +66,15 @@ public static class UInt128Extensions } prime = smallPrimes[i]; - while (order % prime == 0UL) - { - temp = order / prime; - if (temp.PowModWithCycle(q, cycle) == one) - { - order = temp; - } + while (order % prime == 0UL) + { + temp = order / prime; + // TODO: Switch this divisor-order powmod to the ProcessEightBitWindows helper so the + // cycle factoring loop benefits from the faster windowed pow2 ladder measured in CPU benchmarks. + if (temp.PowModWithCycle(q, cycle) == one) + { + order = temp; + } else { break; @@ -103,19 +105,22 @@ public static class UInt128Extensions uint[] smallPrimes = PrimesGenerator.SmallPrimes; ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; ulong i, smallPrimesCount = (ulong)smallPrimes.Length; - for (i = 0UL; i < smallPrimesCount; i++) - { - if (smallPrimesPow2[i] > n) - { - break; - } + for (i = 0UL; i < smallPrimesCount; i++) + { + if (smallPrimesPow2[i] > n) + { + break; + } - p = smallPrimes[i]; - if (n % p == zero) - { - return n == p; - } - } + p = smallPrimes[i]; + // TODO: Replace this direct `%` test with the shared divisor-cycle filter once the + // UInt128 path is wired into the cached cycle tables so wide candidates skip the slow + // modulo checks during primality pre-filtering. + if (n % p == zero) + { + return n == p; + } + } return true; } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index a2eee9b..2bf6fac 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -29,6 +29,9 @@ public static class ULongExtensions } prime = smallPrimes[i]; + // TODO: Replace this `%` driven factor peeling with the divisor-cycle aware + // factoring helper so large orders reuse the cached remainders highlighted in + // the latest divisor-cycle benchmarks instead of recomputing slow modulo checks. while (order % prime == 0UL) { temp = order / prime; @@ -77,6 +80,9 @@ public static class ULongExtensions } p = smallPrimes[i]; + // TODO: Swap this modulo check for the shared small-prime cycle filter once the + // divisor-cycle cache is mandatory, matching the PrimeTester improvements noted in + // the CPU sieve benchmarks. if ((n % p) == 0UL) { return n == p; @@ -86,6 +92,8 @@ public static class ULongExtensions return true; } + // TODO: Swap the `% 5` branches with the Mod5 lookup helper once it is wired into + // production, matching the Mod8/Mod10 benchmark optimizations for hot CLI loops. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod10(this ulong value) => (value & 1UL) == 0UL ? (value % 5UL) switch @@ -108,51 +116,83 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod128(this ulong value) => value & 127UL; + // TODO: Promote the cached Mod5/Mod3 helpers measured fastest in the CLI benchmarks to + // remove the raw modulo operations in this combined remainder computation. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void Mod10_8_5_3(this ulong value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3) { mod8 = value & 7UL; + // TODO: Swap these modulo operations for the dedicated Mod5/Mod3 helpers once the + // benchmarked lookup-based reducers land so this hot helper avoids `%` entirely. + mod5 = value % 5UL; + mod3 = value % 3UL; - byte current; - uint byteSum = 0U; - do - { - current = (byte)value; - byteSum += current; - value >>= 8; - } - while (value != 0UL); - - mod3 = byteSum % 3U; - mod5 = byteSum % 5U; mod10 = (mod8 & 1UL) == 0UL - ? mod5 switch - { - 0UL => 0UL, - 1UL => 6UL, - 2UL => 2UL, - 3UL => 8UL, - _ => 4UL, - } - : mod5 switch - { - 0UL => 5UL, - 1UL => 1UL, - 2UL => 7UL, - 3UL => 3UL, - _ => 9UL, - }; + ? mod5 switch + { + 0UL => 0UL, + 1UL => 6UL, + 2UL => 2UL, + 3UL => 8UL, + _ => 4UL, + } + : mod5 switch + { + 0UL => 5UL, + 1UL => 1UL, + 2UL => 7UL, + 3UL => 3UL, + _ => 9UL, + }; } + // TODO: Inline the Mod5/Mod3 lookup tables here once they are shared so the stepping helper + // mirrors the benchmarked no-modulo variant planned for the CLI sieve. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void Mod10_8_5_3Steps(this ulong value, out ulong step10, out ulong step8, out ulong step5, out ulong step3) { - value.Mod10_8_5_3(out step10, out step8, out step5, out step3); + ulong mod8 = value & 7UL; + // TODO: Route these modulo computations through the cached Mod5/Mod3 tables when the + // lookup implementation from the CLI benchmarks is promoted to production. + ulong mod5 = value % 5UL; + ulong mod3 = value % 3UL; + ulong mod10 = (mod8 & 1UL) == 0UL + ? mod5 switch + { + 0UL => 0UL, + 1UL => 6UL, + 2UL => 2UL, + 3UL => 8UL, + _ => 4UL, + } + : mod5 switch + { + 0UL => 5UL, + 1UL => 1UL, + 2UL => 7UL, + 3UL => 3UL, + _ => 9UL, + }; - step10 = (step10 << 1).Mod10(); - step8 = (step8 << 1) & 7UL; - step5 = (step5 << 1) % 5UL; - step3 = (step3 << 1) % 3UL; + step10 = mod10 + mod10; + if (step10 >= 10UL) + { + step10 -= 10UL; + } + + step8 = (mod8 + mod8) & 7UL; + + step5 = mod5 + mod5; + if (step5 >= 5UL) + { + step5 -= 5UL; + } + + step3 = mod3 + mod3; + if (step3 >= 3UL) + { + step3 -= 3UL; + } } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -161,6 +201,9 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulHigh(this ulong x, ulong y) { + // TODO: Investigate replacing this manual decomposition with the UInt128-based implementation + // for CPU callers; the latest benchmarks show the intrinsic path is an order of magnitude + // faster, while GPU code can keep using GpuUInt128.MulHigh. ulong xLow = (uint)x; ulong xHigh = x >> 32; ulong yLow = (uint)y; @@ -209,11 +252,15 @@ 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) { + // 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); } @@ -221,6 +268,8 @@ public static class ULongExtensions [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); } @@ -259,6 +308,8 @@ public static class ULongExtensions { if ((remainingExponent & 1UL) != 0UL) { + // TODO: Switch to the ProcessEightBitWindows helper once the scalar implementation ships so CPU callers + // inherit the same 2× speedup that the benchmarked windowed ladder delivered for large exponents. result = result.MontgomeryMultiply(baseVal, modulus, nPrime); } @@ -291,6 +342,8 @@ public static class ULongExtensions { if ((remainingExponent & 1UL) != 0UL) { + // TODO: Route this Montgomery-domain variant through the same windowed pow2 helper once available so delta + // stepping matches the ProcessEightBitWindows gains highlighted by GpuPow2ModBenchmarks. result = result.MontgomeryMultiply(baseVal, modulus, divisor.NPrime); } @@ -335,6 +388,8 @@ public static class ULongExtensions UInt128 result = UInt128.One; ulong exponentLoopIndex = 0UL; + // TODO: Port this scalar PowMod fallback to the ProcessEightBitWindows helper so CPU callers get the + // eight-bit window wins measured against the classic square-and-subtract implementation. // Return 1 because 2^0 = 1 if (exponent == 0UL) return result; @@ -385,6 +440,8 @@ public static class ULongExtensions UInt128 result = UInt128.One; ulong exponentLoopIndex = 0UL; + // TODO: Wire this cycle-aware overload into the ProcessEightBitWindows helper so the reduced exponent path + // inherits the faster windowed pow2 routine highlighted in the Pow2Montgomery benchmarks. // Return 1 because 2^0 = 1 if (exponent == 0UL) return result; @@ -436,6 +493,8 @@ public static class ULongExtensions UInt128 result = UInt128.One; ulong exponentLoopIndex = 0UL; + // TODO: Replace this UInt128-cycle overload with the ProcessEightBitWindows helper so large-exponent CPU scans + // reuse the faster windowed pow2 ladder instead of the manual rotation loop measured to lag behind in benchmarks. // Return 1 because 2^0 = 1 if (exponent == UInt128.Zero) return result; @@ -490,6 +549,8 @@ public static class ULongExtensions zero = UInt128.Zero; ulong exponentLoopIndex = 0UL; + // TODO: Migrate this UInt128 exponent overload to ProcessEightBitWindows so the large-cycle reductions drop the + // slow manual loop that underperforms the windowed pow2 helper in the Pow2 benchmark suite. // Return 1 because 2^0 = 1 if (exponent == zero) return result;