Consolidate MulMod helpers across CPU factoring path
This commit is contained in:
Marek Lemańczyk committed
1 parent
bfc6605e97
commit
34240ab3e9
4 files changed
+173
-176
No files matched your search
@@ -105,7 +105,7 @@ public static class MersennePrimeFactorTester
|
||||
}
|
||||
|
||||
candidate = order / prime;
|
||||
if (ModPow128(2, candidate, q, ct) == 1)
|
||||
if (UInt128Numbers.Two.ModPow(candidate, q, ct) == UInt128.One)
|
||||
{
|
||||
order = candidate;
|
||||
}
|
||||
@@ -130,81 +130,6 @@ public static class MersennePrimeFactorTester
|
||||
return order;
|
||||
}
|
||||
|
||||
[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;
|
||||
ulong baseValue = value % modulus;
|
||||
ulong exp = exponent;
|
||||
|
||||
while (exp != 0UL)
|
||||
{
|
||||
if (ct.IsCancellationRequested)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
|
||||
if ((exp & 1UL) != 0UL)
|
||||
{
|
||||
result = MulMod64(result, baseValue, modulus);
|
||||
}
|
||||
|
||||
if (ct.IsCancellationRequested)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
|
||||
baseValue = MulMod64(baseValue, baseValue, modulus);
|
||||
exp >>= 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 +316,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 +459,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)
|
||||
{
|
||||
|
||||
@@ -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;
|
||||
}
|
||||
}
|
||||
@@ -1,5 +1,6 @@
|
||||
using System.Numerics;
|
||||
using System.Runtime.CompilerServices;
|
||||
using System.Threading;
|
||||
using PerfectNumbers.Core.Gpu;
|
||||
|
||||
namespace PerfectNumbers.Core;
|
||||
@@ -559,6 +560,43 @@ public static class UInt128Extensions
|
||||
public static ulong Mod128(this UInt128 value) => (ulong)value & 127UL;
|
||||
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static UInt128 ModPow(this UInt128 value, UInt128 exponent, UInt128 modulus, CancellationToken ct)
|
||||
{
|
||||
UInt128 zero = UInt128.Zero;
|
||||
UInt128 one = UInt128.One;
|
||||
if (modulus <= one)
|
||||
{
|
||||
return zero;
|
||||
}
|
||||
|
||||
UInt128 result = UInt128.One;
|
||||
UInt128 baseValue = value % modulus;
|
||||
|
||||
while (exponent != zero)
|
||||
{
|
||||
if (ct.IsCancellationRequested)
|
||||
{
|
||||
return zero;
|
||||
}
|
||||
|
||||
if ((exponent & one) != zero)
|
||||
{
|
||||
result = result.MulMod(baseValue, modulus);
|
||||
}
|
||||
|
||||
if (ct.IsCancellationRequested)
|
||||
{
|
||||
return zero;
|
||||
}
|
||||
|
||||
baseValue = baseValue.MulMod(baseValue, modulus);
|
||||
exponent >>= 1;
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static UInt128 ModPow(this UInt128 value, UInt128 exponent, UInt128 modulus)
|
||||
{
|
||||
@@ -630,6 +668,46 @@ public static class UInt128Extensions
|
||||
return left;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static UInt128 AddMod(this UInt128 value, UInt128 addend, UInt128 modulus)
|
||||
{
|
||||
if (modulus <= UInt128.One)
|
||||
{
|
||||
return UInt128.Zero;
|
||||
}
|
||||
|
||||
value %= modulus;
|
||||
addend %= modulus;
|
||||
if (addend == UInt128.Zero)
|
||||
{
|
||||
return value;
|
||||
}
|
||||
|
||||
UInt128 threshold = modulus - addend;
|
||||
if (value >= threshold)
|
||||
{
|
||||
return value - threshold;
|
||||
}
|
||||
|
||||
return value + addend;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static UInt128 SubtractOneMod(this UInt128 value, UInt128 modulus)
|
||||
{
|
||||
if (modulus <= UInt128.One)
|
||||
{
|
||||
return UInt128.Zero;
|
||||
}
|
||||
|
||||
if (value == UInt128.Zero)
|
||||
{
|
||||
return modulus - UInt128.One;
|
||||
}
|
||||
|
||||
return value - UInt128.One;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static UInt128 Pow(this UInt128 value, ulong exponent)
|
||||
{
|
||||
|
||||
@@ -455,13 +455,66 @@ public static partial class ULongExtensions
|
||||
return result;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static ulong AddMod64(this ulong value, ulong addend, ulong modulus)
|
||||
{
|
||||
if (modulus <= 1UL)
|
||||
{
|
||||
return 0UL;
|
||||
}
|
||||
|
||||
UInt128 sum = (UInt128)value + addend;
|
||||
UInt128 mod = modulus;
|
||||
if (sum >= mod)
|
||||
{
|
||||
sum -= mod;
|
||||
if (sum >= mod)
|
||||
{
|
||||
sum -= mod;
|
||||
}
|
||||
}
|
||||
|
||||
return (ulong)sum;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private static ulong FoldMontgomery(ulong value, ulong modulus)
|
||||
{
|
||||
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)
|
||||
{
|
||||
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 +530,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)
|
||||
|
||||
Reference in new issue
Block a user