diff --git a/PerfectNumbers.Core/MutableUInt128.cs b/PerfectNumbers.Core/MutableUInt128.cs new file mode 100644 index 0000000..1b53b11 --- /dev/null +++ b/PerfectNumbers.Core/MutableUInt128.cs @@ -0,0 +1,161 @@ +using System; +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +/// +/// Mutable helper for 128-bit arithmetic derived from the .NET runtime UInt128 implementation. +/// Portions adapted from https://github.com/dotnet/runtime/blob/v8.0.0/src/libraries/System.Private.CoreLib/src/System/UInt128.cs. +/// +internal struct MutableUInt128 +{ + public ulong High; + public ulong Low; + + public MutableUInt128(ulong value) + { + High = 0UL; + Low = value; + } + + public MutableUInt128(ulong high, ulong low) + { + High = high; + Low = low; + } + + public MutableUInt128(UInt128 value) + { + High = (ulong)(value >> 64); + Low = (ulong)value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Set(ulong value) + { + High = 0UL; + Low = value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Set(UInt128 value) + { + High = (ulong)(value >> 64); + Low = (ulong)value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Set(ulong high, ulong low) + { + High = high; + Low = low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public UInt128 ToUInt128() => ((UInt128)High << 64) | Low; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Add(ulong value) + { + ulong low = Low + value; + if (low < value) + { + High++; + } + + Low = low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Add(in MutableUInt128 value) + { + ulong low = Low + value.Low; + ulong carry = low < Low ? 1UL : 0UL; + ulong high = High + value.High + carry; + + Low = low; + High = high; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Subtract(ulong value) + { + ulong low = Low; + ulong result = low - value; + ulong borrow = result > low ? 1UL : 0UL; + + Low = result; + High -= borrow; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Subtract(in MutableUInt128 value) + { + ulong low = Low; + ulong result = low - value.Low; + ulong borrow = low < value.Low ? 1UL : 0UL; + + Low = result; + High = High - value.High - borrow; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Multiply(ulong value) + { + ulong currentHigh = High; + ulong currentLow = Low; + + ulong carry = Multiply64(currentLow, value, out ulong newLow); + Multiply64(currentHigh, value, out ulong shiftedHigh); + ulong newHigh = shiftedHigh + carry; + + Low = newLow; + High = newHigh; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Multiply64(ulong left, ulong right, out ulong low) + { + const ulong Mask32 = 0xFFFF_FFFFUL; + + ulong leftLow = left & Mask32; + ulong leftHigh = left >> 32; + ulong rightLow = right & Mask32; + ulong rightHigh = right >> 32; + + ulong lowLow = leftLow * rightLow; + ulong lowHigh = leftLow * rightHigh; + ulong highLow = leftHigh * rightLow; + ulong highHigh = leftHigh * rightHigh; + + ulong carry = (lowLow >> 32) + (lowHigh & Mask32) + (highLow & Mask32); + ulong high = highHigh + (lowHigh >> 32) + (highLow >> 32) + (carry >> 32); + + low = ((carry & Mask32) << 32) | (lowLow & Mask32); + return high; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void MultiplyAdd(ulong multiplier, ulong addend) + { + Multiply(multiplier); + Add(addend); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong Mod(ulong modulus) + { + if (modulus == 0UL) + { + throw new DivideByZeroException(); + } + + if (High == 0UL) + { + return Low % modulus; + } + + return (ulong)(ToUInt128() % modulus); + } +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 1312815..f139cee 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -1,7 +1,6 @@ using System.Buffers; using System.Buffers.Binary; using System.Diagnostics; -using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using PerfectNumbers.Core.Gpu; @@ -73,6 +72,9 @@ internal static partial class PrimeOrderCalculator [ThreadStatic] private static int s_divisorCycleSuppressionDepth; + [ThreadStatic] + private static MutableUInt128 s_pollardRhoPolynomialBuffer; + private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; private readonly struct Pow2ModeScope @@ -314,69 +316,84 @@ internal static partial class PrimeOrderCalculator return order; } - [ThreadStatic] - private static ArrayPool? _factorPool; - private static ArrayPool FactorPool => _factorPool ??= ArrayPool.Create(); + [ThreadStatic] + private static ArrayPool? _factorPool; + private static ArrayPool FactorPool => _factorPool ??= ArrayPool.Create(); - [ThreadStatic] - private static ArrayPool? _ulongPool; - private static ArrayPool ULongPool => _ulongPool ??= ArrayPool.Create(); + private const int FactorEntryStackThreshold = 32; - private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - ArrayPool factorPool = FactorPool; + [ThreadStatic] + private static ArrayPool? _ulongPool; + private static ArrayPool ULongPool => _ulongPool ??= ArrayPool.Create(); - FactorEntry[]? tempArray = null; - try - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - tempArray = factorPool.Rent(length + 1); - Span buffer = tempArray; - factorSpan.CopyTo(buffer); + private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + ArrayPool factorPool = FactorPool; - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None); - if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime) - { - buffer[length] = new FactorEntry(factors.Cofactor); - length++; - } + FactorEntry[]? tempArray = null; + try + { + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + int capacity = length + 1; - buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value)); + Span buffer; + if (capacity <= FactorEntryStackThreshold) + { + buffer = stackalloc FactorEntry[capacity]; + } + else + { + tempArray = factorPool.Rent(capacity); + buffer = tempArray.AsSpan(0, capacity); + } - for (int i = 0; i < length; i++) - { - ulong primeFactor = buffer[i].Value; - int exponent = buffer[i].Exponent; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } + factorSpan.CopyTo(buffer); - ulong reduced = order / primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - order = reduced; - continue; - } + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None); + if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime) + { + buffer[length] = new FactorEntry(factors.Cofactor); + length++; + } - break; - } - } + buffer = buffer.Slice(0, length); // Reuse the same span restricted to populated entries. + buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - return order; - } - finally - { - if (tempArray is not null) - { - factorPool.Return(tempArray, clearArray: false); - } - } - } + for (int i = 0; i < buffer.Length; i++) + { + FactorEntry entry = buffer[i]; + ulong primeFactor = entry.Value; + int exponent = entry.Exponent; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong reduced = order / primeFactor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + order = reduced; + continue; + } + + break; + } + } + + return order; + } + finally + { + if (tempArray is not null) + { + factorPool.Return(tempArray, clearArray: false); + } + } + } private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, ref PartialFactorResult orderFactors) { @@ -1311,8 +1328,11 @@ internal static partial class PrimeOrderCalculator [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) { - UInt128 value = (UInt128)x * x + c; - return (ulong)(value % modulus); + ref MutableUInt128 polynomial = ref s_pollardRhoPolynomialBuffer; + polynomial.Set(x); + polynomial.Multiply(x); + polynomial.Add(c); + return polynomial.Mod(modulus); } private static ulong BinaryGcd(ulong a, ulong b) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 8c08577..a2ed557 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -346,48 +346,67 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount); + + Span oddPowerStorage; + ulong[]? pooledOddPowers = null; + + if (oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount) + { + oddPowerStorage = stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]; + } + else if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold) + { + oddPowerStorage = new ulong[oddPowerCount]; + } + else + { + pooledOddPowers = ArrayPool.Shared.Rent(oddPowerCount); + oddPowerStorage = pooledOddPowers; + } + + Span oddPowers = oddPowerStorage[..oddPowerCount]; + InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers); int index = bitLength - 1; - while (index >= 0) - { - if (((exponent >> index) & 1UL) == 0UL) - { - result = result.MontgomeryMultiply(result, modulus, nPrime); - index--; - continue; - } + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + index--; + continue; + } - int windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } - while (((exponent >> windowStart) & 1UL) == 0UL) - { - windowStart++; - } + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } - int windowLength = index - windowStart + 1; - for (int square = 0; square < windowLength; square++) - { - result = result.MontgomeryMultiply(result, modulus, nPrime); - } + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + } - ulong mask = (1UL << windowLength) - 1UL; - ulong windowValue = (exponent >> windowStart) & mask; - int tableIndex = (int)((windowValue - 1UL) >> 1); - ulong multiplier = oddPowers[tableIndex]; - result = result.MontgomeryMultiply(multiplier, modulus, nPrime); + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MontgomeryMultiply(multiplier, modulus, nPrime); - index = windowStart - 1; - } - - if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold) - { - ArrayPool.Shared.Return(oddPowers); - } + index = windowStart - 1; + } + + if (pooledOddPowers is not null) + { + ArrayPool.Shared.Return(pooledOddPowers); + } if (keepMontgomery) { @@ -533,29 +552,25 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span destination) { - ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold - ? new ulong[oddPowerCount] - : ArrayPool.Shared.Rent(oddPowerCount); - + int oddPowerCount = destination.Length; if (oddPowerCount == 0) { - return oddPowers; + return; } - Span destination = oddPowers.AsSpan(0, oddPowerCount); bool computedOnGpu = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount && MontgomeryOddPowerGpu.TryCompute(divisor, oddPowerCount, destination); if (computedOnGpu) { - return oddPowers; + return; } destination[0] = divisor.MontgomeryTwo; if (oddPowerCount == 1) { - return oddPowers; + return; } ulong square = divisor.MontgomeryTwoSquared; @@ -565,8 +580,6 @@ public static class ULongExtensions current = current.MontgomeryMultiply(square, modulus, nPrime); destination[i] = current; } - - return oddPowers; } [MethodImpl(MethodImplOptions.AggressiveInlining)]