mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #211 from mlemanczyk/codex/analyze-mulmod-calculations-in-evenperfectbitscanner

Browse Source 
Consolidate MulMod helpers across CPU factoring path
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-10-20 11:53:00 +02:00
2 parents 4e553c951b ea2476c749
commit 9c08cf48a7
4 files changed
+172 -157
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/MersennePrimeFactorTester.cs
+24 -70
View File
@@ -105,7 +105,7 @@ public static class MersennePrimeFactorTester
}
candidate = order / prime;
if (ModPow128(2, candidate, q, ct) == 1)
if (ModPowWithCancellation(UInt128Numbers.Two, candidate, q, ct) == UInt128.One)
{
order = candidate;
}
@@ -131,80 +131,42 @@ 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)
private static UInt128 ModPowWithCancellation(UInt128 baseValue, UInt128 exponent, UInt128 modulus, CancellationToken ct)
{
ulong result = 1UL;
ulong baseValue = value % modulus;
ulong exp = exponent;
// Keep this helper local to the factoring path so we can preserve the token-aware early outs
// without carrying a separate overload on the shared UInt128 extensions.
if (modulus <= UInt128.One)
{
return UInt128.Zero;
}
while (exp != 0UL)
UInt128 result = UInt128.One;
UInt128 baseResidue = baseValue % modulus;
while (exponent != UInt128.Zero)
{
if (ct.IsCancellationRequested)
{
return 0UL;
return UInt128.Zero;
}
if ((exp & 1UL) != 0UL)
if ((exponent & UInt128.One) != UInt128.Zero)
{
result = MulMod64(result, baseValue, modulus);
result = result.MulMod(baseResidue, modulus);
}
if (ct.IsCancellationRequested)
{
return 0UL;
return UInt128.Zero;
}
baseValue = MulMod64(baseValue, baseValue, modulus);
exp >>= 1;
baseResidue = baseResidue.MulMod(baseResidue, modulus);
exponent >>= 1;
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
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;
UInt128 baseValue = value % modulus;
UInt128 exp = exponent;
while (exp != 0)
{
if (ct.IsCancellationRequested)
{
return 0;
}
if ((exp & 1) != 0)
{
result = MulMod128(result, baseValue, modulus);
}
if (ct.IsCancellationRequested)
{
return 0;
}
baseValue = MulMod128(baseValue, baseValue, modulus);
exp >>= 1;
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 mod)
{
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;
}
private static readonly ConcurrentBag<List<ulong>> ListPool = new();
private static readonly ConcurrentBag<Dictionary<ulong, byte>> DictPool = new();
private static readonly Dictionary<ulong, ulong[]> FactorCache = new();
@@ -391,24 +353,19 @@ public static class MersennePrimeFactorTester
return 0UL;
}
x = (MulMod64(x, x, n) + c) % n; // TODO: Swap the `% n` with the Montgomery folding
// helper from the PollardRho benchmarks once the
// specialized reducer lands.
x = x.MulMod64(x, n).AddMod64(c, n);
if (ct.IsCancellationRequested)
{
return 0UL;
}
y = (MulMod64(y, y, n) + c) % n; // TODO: Use the same Montgomery folding helper here
// to keep the tortoise sequence on the optimized
// path.
y = y.MulMod64(y, n).AddMod64(c, n);
if (ct.IsCancellationRequested)
{
return 0UL;
}
y = (MulMod64(y, y, n) + c) % n; // TODO: Replace this modulo with the optimized helper
// so both steps share the fast reduction.
y = y.MulMod64(y, n).AddMod64(c, n);
diff = x > y ? x - y : y - x;
if (ct.IsCancellationRequested)
{
@@ -539,22 +496,19 @@ public static class MersennePrimeFactorTester
return 0;
}
x = (MulMod128(x, x, n) + c) % n; // TODO: Swap this modulo for the UInt128 Montgomery
// reducer once the optimized helper lands.
x = x.MulMod(x, n).AddMod(c, n);
if (ct.IsCancellationRequested)
{
return 0;
}
y = (MulMod128(y, y, n) + c) % n; // TODO: Use the same optimized reducer for the
// tortoise step.
y = y.MulMod(y, n).AddMod(c, n);
if (ct.IsCancellationRequested)
{
return 0;
}
y = (MulMod128(y, y, n) + c) % n; // TODO: Replace with the optimized reducer so both
// steps avoid BigInteger-based modulo.
y = y.MulMod(y, n).AddMod(c, n);
diff = x > y ? x - y : y - x;
if (ct.IsCancellationRequested)
{
PerfectNumbers.Core/ModResidueTracker.cs
+26 -80
View File
@@ -44,7 +44,7 @@ public sealed class ModResidueTracker
if (divisor == UInt128.Zero)
{
throw new ArgumentOutOfRangeException(nameof(divisor));
}
}
EnsureAt(number);
@@ -58,15 +58,15 @@ public sealed class ModResidueTracker
if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended)
{
_lastAppended = divisor;
}
}
else if (!_hasLastAppended)
{
_hasLastAppended = true;
_lastAppended = divisor;
}
}
return residue == UInt128.Zero;
}
}
return _residues[idx] == UInt128.Zero;
}
@@ -90,7 +90,7 @@ public sealed class ModResidueTracker
EnsureAt(number);
// After EnsureAt we can keep _cursor (still valid monotonic pointer)
if (_cursor > divisors.Count) _cursor = divisors.Count;
}
}
UInt128 residue;
// Fast path: append to tail when candidates are strictly increasing.
@@ -105,21 +105,21 @@ public sealed class ModResidueTracker
_cursor = divisors.Count; // position after last
divisible = residue == UInt128.Zero;
return true;
}
}
// Merge pointer: advance until current >= divisor
var cursor = _cursor;
while (cursor < divisors.Count && divisors[cursor] < divisor)
{
cursor++;
}
}
if (cursor < divisors.Count && divisors[cursor] == divisor)
{
divisible = residues[cursor] == UInt128.Zero;
_cursor = cursor + 1; // position after matched
return true;
}
}
// Insert at cursor to keep order
@@ -134,7 +134,7 @@ public sealed class ModResidueTracker
{
_hasLastAppended = true;
_lastAppended = divisor;
}
}
return true;
}
@@ -150,17 +150,17 @@ public sealed class ModResidueTracker
// Initialize residues lazily upon first check for each divisor.
_initialized = true;
return;
}
}
if (number == _currentNumber)
{
return;
}
}
if (number < _currentNumber)
{
throw new InvalidOperationException("Number must be non-decreasing.");
}
}
UInt128 delta = number - _currentNumber;
List<UInt128> divisors = _divisors;
@@ -169,7 +169,7 @@ public sealed class ModResidueTracker
{
_currentNumber = number;
return;
}
}
// Advance all tracked residues from _currentNumber to 'number'.
// For Identity: r' = (r + delta) mod d
@@ -206,17 +206,15 @@ public sealed class ModResidueTracker
// 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);
UInt128 powDelta = UInt128Numbers.Two.ModPow(delta, d);
// r' = pow2Delta*r + (pow2Delta - 1) mod d
value = MulMod128(value, residues[i], d) + value - UInt128.One;
if (value >= d)
{
value -= d;
}
value = powDelta.MulMod(residues[i], d);
value = value.AddMod(powDelta, d);
value = value.SubtractOneMod(d);
residues[i] = value;
}
}
residues[i] = value;
}
}
_currentNumber = number;
}
@@ -231,12 +229,12 @@ public sealed class ModResidueTracker
if (list[mid] < value)
{
lo = mid + 1;
}
}
else
{
hi = mid;
}
}
}
}
return lo;
}
@@ -264,70 +262,18 @@ public sealed class ModResidueTracker
// instead of recomputing powmods for every divisor; the divisor-cycle benchmarks showed large wins once the
// cached orders were used across scans, and a miss should trigger the configured device to compute the cycle
// immediately without storing it back or requesting additional cache blocks.
UInt128 pow = PowMod128(UInt128Numbers.Two, number, divisor);
UInt128 pow = UInt128Numbers.Two.ModPow(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;
ulong ord = divisor.CalculateOrder();
if (ord == 0UL || divisor % ord != UInt128.Zero)
if (ord == 0UL || number % (UInt128)ord != UInt128.Zero)
return UInt128.One; // does NOT divide
return UInt128.Zero; // q really divides M_p
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
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;
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)
{
value = MulMod128(value, value, modulus);
}
}
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.
// 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)
{
if ((b & UInt128.One) != UInt128.Zero)
{
result += result < modulus ? a : a - modulus;
}
a = a < modulus ? a << 1 : (a << 1) - modulus;
b >>= 1;
}
return result;
}
}
PerfectNumbers.Core/UInt128Extensions.cs
+52
View File
@@ -1,5 +1,6 @@
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Threading;
using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core;
@@ -630,6 +631,57 @@ public static class UInt128Extensions
return left;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 AddMod(this UInt128 value, UInt128 addend, UInt128 modulus)
{
// All current CPU callers (Pollard Rho and the Mersenne residue tracker) only supply moduli
// of at least three, so leave this guard documented but disabled.
// if (modulus <= UInt128.One)
// {
// return UInt128.Zero;
// }
// Both operands arrive pre-reduced on these paths, letting us skip the redundant folds the old
// helper performed.
// value %= modulus;
// addend %= modulus;
// Pollard Rho keeps its polynomial constant inside [1, modulus - 1], and the Mersenne pow-delta
// term never hits zero for odd moduli, leaving this branch inactive outside targeted tests.
// if (addend == UInt128.Zero)
// {
// return value;
// }
UInt128 threshold = modulus - addend;
if (value >= threshold)
{
return value - threshold;
}
UInt128 sum = value + addend;
// The sum stays below the modulus under the threshold guard above, so no wraparound occurs.
return sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 SubtractOneMod(this UInt128 value, UInt128 modulus)
{
// Only the Mersenne residue tracker calls this helper, and it never observes moduli below three,
// so leave this guard documented but inactive.
// if (modulus <= UInt128.One)
// {
// return UInt128.Zero;
// }
if (value == UInt128.Zero)
{
// Wrapping occurs when the prior AddMod call produced one, so step back to modulus - 1.
return modulus - UInt128.One;
}
return value - UInt128.One;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow(this UInt128 value, ulong exponent)
{
PerfectNumbers.Core/ULongExtensions.cs
+70 -7
View File
@@ -455,13 +455,73 @@ public static partial class ULongExtensions
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong AddMod64(this ulong value, ulong addend, ulong modulus)
{
// Every CPU caller (Pollard Rho and the ModPow64 ladder) supplies moduli of at least three, so leave this guard
// documented but disabled.
// if (modulus <= 1UL)
// {
// return 0UL;
// }
// Both flows keep their operands below the modulus ahead of time, letting us skip the re-fold the old helper performed.
UInt128 sum = (UInt128)value + addend;
if (sum < modulus)
{
return (ulong)sum;
}
sum -= modulus;
// After one subtraction the sum always falls below the modulus on this path, making the old guard redundant.
// if (sum >= modulus)
// {
// sum -= modulus;
// }
return (ulong)sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong FoldMontgomery(ulong value, ulong modulus)
{
// Pollard Rho and ModPow64 both feed moduli of at least three, so retain the documentation and skip the runtime guard.
// if (modulus <= 1UL)
// {
// return 0UL;
// }
if (value < modulus)
{
return value;
}
ulong folded = 0UL;
int leadingZeros = BitOperations.LeadingZeroCount(value);
int bitIndex = 63 - leadingZeros;
for (; bitIndex >= 0; bitIndex--)
{
folded = folded.AddMod64(folded, modulus);
if (((value >> bitIndex) & 1UL) != 0UL)
{
folded = folded.AddMod64(1UL, modulus);
}
}
return folded;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus)
{
// Pollard Rho and the Mersenne testers never pass moduli below two, so keep the guard commented for clarity.
// if (modulus <= 1UL)
// {
// return 0UL;
// }
ulong result = 1UL;
// TODO: Replace this `%` with the Montgomery folding helper highlighted in MulMod64Benchmarks so the
// modular exponentiation avoids the slow integer division before the ladder even starts.
value %= modulus;
value = FoldMontgomery(value, modulus);
while (exponent != 0UL)
{
@@ -477,10 +537,13 @@ public static partial class ULongExtensions
return result;
}
[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)]
// Reuse the direct UInt128 reduction path that dominated MulMod64Benchmarks to avoid
// the redundant operand modulo operations in the older fallback.
public static ulong MulMod64(this ulong a, ulong b, ulong modulus)
{
return (ulong)(((UInt128)a * b) % modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)