diff --git a/PerfectNumbers.Core/Cpu/MersennePollardBrent.cs b/PerfectNumbers.Core/Cpu/MersennePollardBrent.cs index d8ffba2..7e0c587 100644 --- a/PerfectNumbers.Core/Cpu/MersennePollardBrent.cs +++ b/PerfectNumbers.Core/Cpu/MersennePollardBrent.cs @@ -2,17 +2,32 @@ using System; using System.Buffers; using System.Numerics; using System.Security.Cryptography; +using PerfectNumbers.Core; +using UInt128 = System.UInt128; namespace PerfectNumbers.Core.Cpu; internal sealed class MersennePollardBrent { + private enum ArithmeticMode + { + Unsupported, + UInt64, + UInt128, + BigInteger, + } + + private readonly ArithmeticMode _mode; private readonly int _bitLength; private readonly int _byteLength; private readonly int _bitRemainder; private readonly byte _maskByte; private readonly BigInteger _modulus; private readonly BigInteger _mask; + private readonly ulong _modulus64; + private readonly ulong _mask64; + private readonly UInt128 _modulus128; + private readonly UInt128 _mask128; private readonly int _batchSize; private readonly bool _isSupported; @@ -20,6 +35,7 @@ internal sealed class MersennePollardBrent { if (exponent < 3UL || exponent > int.MaxValue) { + _mode = ArithmeticMode.Unsupported; _isSupported = false; _bitLength = 0; _byteLength = 0; @@ -27,6 +43,10 @@ internal sealed class MersennePollardBrent _maskByte = 0; _modulus = BigInteger.Zero; _mask = BigInteger.Zero; + _modulus64 = 0UL; + _mask64 = 0UL; + _modulus128 = UInt128.Zero; + _mask128 = UInt128.Zero; _batchSize = 0; return; } @@ -36,9 +56,57 @@ internal sealed class MersennePollardBrent _byteLength = (_bitLength + 7) / 8; _bitRemainder = _bitLength & 7; _maskByte = _bitRemainder == 0 ? (byte)0xFF : (byte)((1 << _bitRemainder) - 1); + _batchSize = Math.Clamp(batchSize, 8, 1024); + + if (_bitLength <= 64) + { + _mode = ArithmeticMode.UInt64; + if (_bitLength == 64) + { + _modulus64 = ulong.MaxValue; + _mask64 = ulong.MaxValue; + } + else + { + _modulus64 = (1UL << _bitLength) - 1UL; + _mask64 = _modulus64; + } + + _modulus128 = UInt128.Zero; + _mask128 = UInt128.Zero; + _modulus = BigInteger.Zero; + _mask = BigInteger.Zero; + return; + } + + if (_bitLength <= 128) + { + _mode = ArithmeticMode.UInt128; + if (_bitLength == 128) + { + _modulus128 = UInt128.MaxValue; + _mask128 = UInt128.MaxValue; + } + else + { + _modulus128 = (UInt128.One << _bitLength) - UInt128.One; + _mask128 = _modulus128; + } + + _modulus64 = 0UL; + _mask64 = 0UL; + _modulus = BigInteger.Zero; + _mask = BigInteger.Zero; + return; + } + + _mode = ArithmeticMode.BigInteger; _modulus = (BigInteger.One << _bitLength) - BigInteger.One; _mask = _modulus; - _batchSize = Math.Clamp(batchSize, 8, 1024); + _modulus64 = 0UL; + _mask64 = 0UL; + _modulus128 = UInt128.Zero; + _mask128 = UInt128.Zero; } public bool IsSupported => _isSupported; @@ -63,14 +131,222 @@ internal sealed class MersennePollardBrent return false; } + return _mode switch + { + ArithmeticMode.UInt64 => TryFindFactorUInt64(allowedMax, maxRestarts, maxIterationsPerRestart, out factor, out iterationCount), + ArithmeticMode.UInt128 => TryFindFactorUInt128(allowedMax, maxRestarts, maxIterationsPerRestart, out factor, out iterationCount), + ArithmeticMode.BigInteger => TryFindFactorBigInteger(allowedMax, maxRestarts, maxIterationsPerRestart, out factor, out iterationCount), + _ => false, + }; + } + + private bool TryFindFactorUInt64( + ulong allowedMax, + int maxRestarts, + int maxIterationsPerRestart, + out ulong factor, + out ulong iterationCount) + { + factor = 0UL; + iterationCount = 0UL; + int restarts = Math.Max(1, maxRestarts); int iterationBudget = Math.Max(_batchSize, maxIterationsPerRestart); for (int restart = 0; restart < restarts; restart++) { - BigInteger x = NextRandomSeed(); + ulong x = NextRandomSeed64(); + ulong y = x; + ulong c = NextRandomConstant64(); + ulong q = 1UL; + ulong g = 1UL; + ulong ys = 0UL; + int r = 1; + ulong localIterations = 0UL; + + while (g == 1UL && localIterations < (ulong)iterationBudget) + { + x = y; + for (int i = 0; i < r && localIterations < (ulong)iterationBudget; i++) + { + y = Iterate64(y, c, ref localIterations); + } + + if (localIterations >= (ulong)iterationBudget) + { + break; + } + + int k = 0; + while (k < r && g == 1UL && localIterations < (ulong)iterationBudget) + { + ys = y; + int t = Math.Min(_batchSize, r - k); + for (int i = 0; i < t; i++) + { + y = Iterate64(y, c, ref localIterations); + ulong diff = AbsSub64(y, x); + if (diff == 0UL) + { + continue; + } + + q = MulMod64(q, diff); + } + + g = BinaryGcd64(q, _modulus64); + k += t; + } + + r <<= 1; + } + + if (g == _modulus64) + { + do + { + ys = Iterate64(ys, c, ref localIterations); + ulong diff = AbsSub64(ys, x); + if (diff == 0UL) + { + continue; + } + + g = BinaryGcd64(diff, _modulus64); + } + while (g == 1UL); + } + + iterationCount += localIterations; + + if (g == 1UL || g == _modulus64) + { + continue; + } + + if (g > 1UL && g <= allowedMax) + { + factor = g; + return true; + } + } + + return false; + } + + private bool TryFindFactorUInt128( + ulong allowedMax, + int maxRestarts, + int maxIterationsPerRestart, + out ulong factor, + out ulong iterationCount) + { + factor = 0UL; + iterationCount = 0UL; + + int restarts = Math.Max(1, maxRestarts); + int iterationBudget = Math.Max(_batchSize, maxIterationsPerRestart); + UInt128 allowedMaxWide = (UInt128)allowedMax; + + for (int restart = 0; restart < restarts; restart++) + { + UInt128 x = NextRandomSeed128(); + UInt128 y = x; + UInt128 c = NextRandomConstant128(); + UInt128 q = UInt128.One; + UInt128 g = UInt128.One; + UInt128 ys = UInt128.Zero; + int r = 1; + ulong localIterations = 0UL; + + while (g == UInt128.One && localIterations < (ulong)iterationBudget) + { + x = y; + for (int i = 0; i < r && localIterations < (ulong)iterationBudget; i++) + { + y = Iterate128(y, c, ref localIterations); + } + + if (localIterations >= (ulong)iterationBudget) + { + break; + } + + int k = 0; + while (k < r && g == UInt128.One && localIterations < (ulong)iterationBudget) + { + ys = y; + int t = Math.Min(_batchSize, r - k); + for (int i = 0; i < t; i++) + { + y = Iterate128(y, c, ref localIterations); + UInt128 diff = AbsSub128(y, x); + if (diff == UInt128.Zero) + { + continue; + } + + q = MulMod128(q, diff); + } + + g = q.BinaryGcd(_modulus128); + k += t; + } + + r <<= 1; + } + + if (g == _modulus128) + { + do + { + ys = Iterate128(ys, c, ref localIterations); + UInt128 diff = AbsSub128(ys, x); + if (diff == UInt128.Zero) + { + continue; + } + + g = diff.BinaryGcd(_modulus128); + } + while (g == UInt128.One); + } + + iterationCount += localIterations; + + if (g == UInt128.One || g == _modulus128) + { + continue; + } + + if (g > UInt128.One && g <= allowedMaxWide) + { + factor = (ulong)g; + return true; + } + } + + return false; + } + + private bool TryFindFactorBigInteger( + ulong allowedMax, + int maxRestarts, + int maxIterationsPerRestart, + out ulong factor, + out ulong iterationCount) + { + factor = 0UL; + iterationCount = 0UL; + + int restarts = Math.Max(1, maxRestarts); + int iterationBudget = Math.Max(_batchSize, maxIterationsPerRestart); + + for (int restart = 0; restart < restarts; restart++) + { + BigInteger x = NextRandomSeedBigInteger(); BigInteger y = x; - BigInteger c = NextRandomConstant(); + BigInteger c = NextRandomConstantBigInteger(); BigInteger q = BigInteger.One; BigInteger g = BigInteger.One; BigInteger ys = BigInteger.Zero; @@ -82,7 +358,7 @@ internal sealed class MersennePollardBrent x = y; for (int i = 0; i < r && localIterations < (ulong)iterationBudget; i++) { - y = Iterate(y, c, ref localIterations); + y = IterateBigInteger(y, c, ref localIterations); } if (localIterations >= (ulong)iterationBudget) @@ -97,14 +373,14 @@ internal sealed class MersennePollardBrent int t = Math.Min(_batchSize, r - k); for (int i = 0; i < t; i++) { - y = Iterate(y, c, ref localIterations); - BigInteger diff = AbsSub(y, x); + y = IterateBigInteger(y, c, ref localIterations); + BigInteger diff = AbsSubBigInteger(y, x); if (diff.IsZero) { continue; } - q = MulMod(q, diff); + q = MulModBigInteger(q, diff); } g = BigInteger.GreatestCommonDivisor(q, _modulus); @@ -118,8 +394,8 @@ internal sealed class MersennePollardBrent { do { - ys = Iterate(ys, c, ref localIterations); - BigInteger diff = AbsSub(ys, x); + ys = IterateBigInteger(ys, c, ref localIterations); + BigInteger diff = AbsSubBigInteger(ys, x); if (diff.IsZero) { continue; @@ -151,25 +427,86 @@ internal sealed class MersennePollardBrent return false; } - private BigInteger Iterate(BigInteger value, BigInteger c, ref ulong iterationCounter) + private ulong Iterate64(ulong value, ulong c, ref ulong iterationCounter) { iterationCounter++; - value = SquareMod(value); - value = AddMod(value, c); + value = SquareMod64(value); + value = AddMod64(value, c); return value; } - private BigInteger AddMod(BigInteger a, BigInteger b) + private UInt128 Iterate128(UInt128 value, UInt128 c, ref ulong iterationCounter) { - return Reduce(a + b); + iterationCounter++; + value = SquareMod128(value); + value = AddMod128(value, c); + return value; } - private BigInteger SquareMod(BigInteger value) + private BigInteger IterateBigInteger(BigInteger value, BigInteger c, ref ulong iterationCounter) { - return Reduce(value * value); + iterationCounter++; + value = SquareModBigInteger(value); + value = AddModBigInteger(value, c); + return value; } - private BigInteger MulMod(BigInteger a, BigInteger b) + private ulong AddMod64(ulong a, ulong b) + { + UInt128 sum = (UInt128)a + b; + return Reduce64(sum); + } + + private UInt128 AddMod128(UInt128 a, UInt128 b) + { + UInt128 sum = a + b; + bool carryOut = sum < a; + if (!carryOut && sum < _modulus128) + { + return sum; + } + + ulong p3 = 0UL; + ulong p2 = carryOut ? 1UL : 0UL; + ulong p1 = (ulong)(sum >> 64); + ulong p0 = (ulong)sum; + return ReduceProductBitwise(p3, p2, p1, p0, _modulus128); + } + + private BigInteger AddModBigInteger(BigInteger a, BigInteger b) + { + return ReduceBigInteger(a + b); + } + + private ulong SquareMod64(ulong value) + { + UInt128 square = (UInt128)value * value; + return Reduce64(square); + } + + private UInt128 SquareMod128(UInt128 value) + { + return MulMod128(value, value); + } + + private BigInteger SquareModBigInteger(BigInteger value) + { + return ReduceBigInteger(value * value); + } + + private ulong MulMod64(ulong a, ulong b) + { + UInt128 product = (UInt128)a * b; + return Reduce64(product); + } + + private UInt128 MulMod128(UInt128 a, UInt128 b) + { + MultiplyFull(a, b, out ulong p3, out ulong p2, out ulong p1, out ulong p0); + return ReduceProductBitwise(p3, p2, p1, p0, _modulus128); + } + + private BigInteger MulModBigInteger(BigInteger a, BigInteger b) { if (b == _modulus) { @@ -181,10 +518,58 @@ internal sealed class MersennePollardBrent a = BigInteger.Zero; } - return Reduce(a * b); + return ReduceBigInteger(a * b); } - private BigInteger Reduce(BigInteger value) + private ulong Reduce64(UInt128 value) + { + if (_bitLength == 64) + { + UInt128 combined = (UInt128)(ulong)value + (value >> 64); + UInt128 modulus = _modulus64; + if (combined >= modulus) + { + combined -= modulus; + if (combined >= modulus) + { + combined -= modulus; + } + } + + return (ulong)combined; + } + + UInt128 mask = _mask64; + UInt128 modulusValue = _modulus64; + UInt128 working = value; + while (true) + { + UInt128 high = working >> _bitLength; + if (high == UInt128.Zero) + { + break; + } + + working = (working & mask) + high; + } + + if (working >= modulusValue) + { + working -= modulusValue; + if (working >= modulusValue) + { + working %= modulusValue; + if (working >= modulusValue) + { + working -= modulusValue; + } + } + } + + return (ulong)working; + } + + private BigInteger ReduceBigInteger(BigInteger value) { if (value.Sign < 0) { @@ -225,31 +610,127 @@ internal sealed class MersennePollardBrent return value; } - private BigInteger NextRandomSeed() + private ulong NextRandomSeed64() + { + ulong value; + do + { + value = NextRandomBelowModulus64(); + } + while (value <= 1UL); + + return value; + } + + private UInt128 NextRandomSeed128() + { + UInt128 value; + do + { + value = NextRandomBelowModulus128(); + } + while (value <= UInt128.One); + + return value; + } + + private BigInteger NextRandomSeedBigInteger() { BigInteger value; do { - value = NextRandomBelowModulus(); + value = NextRandomBelowModulusBigInteger(); } while (value <= BigInteger.One); return value; } - private BigInteger NextRandomConstant() + private ulong NextRandomConstant64() + { + ulong value; + do + { + value = NextRandomBelowModulus64(); + } + while (value == 0UL); + + return value; + } + + private UInt128 NextRandomConstant128() + { + UInt128 value; + do + { + value = NextRandomBelowModulus128(); + } + while (value == UInt128.Zero); + + return value; + } + + private BigInteger NextRandomConstantBigInteger() { BigInteger value; do { - value = NextRandomBelowModulus(); + value = NextRandomBelowModulusBigInteger(); } while (value.IsZero); return value; } - private BigInteger NextRandomBelowModulus() + private ulong NextRandomBelowModulus64() + { + if (_byteLength == 0) + { + return 0UL; + } + + Span buffer = stackalloc byte[_byteLength]; + ulong candidate; + do + { + RandomNumberGenerator.Fill(buffer); + if (_bitRemainder != 0) + { + buffer[_byteLength - 1] &= _maskByte; + } + + candidate = (ulong)ReadUInt128LittleEndian(buffer); + } + while (candidate >= _modulus64); + + return candidate; + } + + private UInt128 NextRandomBelowModulus128() + { + if (_byteLength == 0) + { + return UInt128.Zero; + } + + Span buffer = stackalloc byte[_byteLength]; + UInt128 candidate; + do + { + RandomNumberGenerator.Fill(buffer); + if (_bitRemainder != 0) + { + buffer[_byteLength - 1] &= _maskByte; + } + + candidate = ReadUInt128LittleEndian(buffer); + } + while (candidate >= _modulus128); + + return candidate; + } + + private BigInteger NextRandomBelowModulusBigInteger() { if (_byteLength == 0) { @@ -280,9 +761,166 @@ internal sealed class MersennePollardBrent } } - private static BigInteger AbsSub(BigInteger a, BigInteger b) + private static ulong AbsSub64(ulong a, ulong b) + { + return a >= b ? a - b : b - a; + } + + private static UInt128 AbsSub128(UInt128 a, UInt128 b) + { + return a >= b ? a - b : b - a; + } + + private static BigInteger AbsSubBigInteger(BigInteger a, BigInteger b) { BigInteger diff = a - b; return diff.Sign < 0 ? -diff : diff; } + + private static ulong BinaryGcd64(ulong u, ulong v) + { + if (u == 0UL) + { + return v; + } + + if (v == 0UL) + { + return u; + } + + int shift = BitOperations.TrailingZeroCount(u | v); + u >>= BitOperations.TrailingZeroCount(u); + do + { + v >>= BitOperations.TrailingZeroCount(v); + if (u > v) + { + (u, v) = (v, u); + } + + v -= u; + } + while (v != 0UL); + + return u << shift; + } + + private static void MultiplyFull(UInt128 left, UInt128 right, out ulong p3, out ulong p2, out ulong p1, out ulong p0) + { + ulong leftLow = (ulong)left; + ulong leftHigh = (ulong)(left >> 64); + ulong rightLow = (ulong)right; + ulong rightHigh = (ulong)(right >> 64); + + (ulong High, ulong Low) partial0 = Multiply64(leftLow, rightLow); + (ulong High, ulong Low) partial1 = Multiply64(leftLow, rightHigh); + (ulong High, ulong Low) partial2 = Multiply64(leftHigh, rightLow); + (ulong High, ulong Low) partial3 = Multiply64(leftHigh, rightHigh); + + p0 = partial0.Low; + ulong carry = partial0.High; + ulong sum = carry + partial1.Low; + ulong carryMid = sum < partial1.Low ? 1UL : 0UL; + sum += partial2.Low; + if (sum < partial2.Low) + { + carryMid++; + } + + p1 = sum; + sum = partial1.High + partial2.High; + ulong carryHigh = sum < partial2.High ? 1UL : 0UL; + sum += partial3.Low; + if (sum < partial3.Low) + { + carryHigh++; + } + + sum += carryMid; + if (sum < carryMid) + { + carryHigh++; + } + + p2 = sum; + p3 = partial3.High + carryHigh; + } + + private static (ulong High, ulong Low) Multiply64(ulong left, ulong right) + { + ulong a0 = (uint)left; + ulong a1 = left >> 32; + ulong b0 = (uint)right; + ulong b1 = right >> 32; + + ulong lo = a0 * b0; + ulong mid1 = a1 * b0; + b0 = a0 * b1; + b1 *= a1; + + a0 = (lo >> 32) + (uint)mid1 + (uint)b0; + a1 = (lo & 0xFFFFFFFFUL) | (a0 << 32); + b1 += (mid1 >> 32) + (b0 >> 32) + (a0 >> 32); + + return (b1, a1); + } + + private static UInt128 ReduceProductBitwise(ulong p3, ulong p2, ulong p1, ulong p0, UInt128 modulus) + { + UInt128 remainder = UInt128.Zero; + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p3 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p2 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p1 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p0 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + return remainder; + } + + private static UInt128 ReadUInt128LittleEndian(ReadOnlySpan buffer) + { + UInt128 value = UInt128.Zero; + for (int i = buffer.Length - 1; i >= 0; i--) + { + value <<= 8; + value |= buffer[i]; + } + + return value; + } }