From 2da549fe8cbe3fd5b4c3b777f0b78ce1cde705a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 11 Oct 2025 22:45:46 +0200 Subject: [PATCH 1/3] Reduce GC in CPU odd power ladder --- PerfectNumbers.Core/ULongExtensions.cs | 105 ++++++++++++++----------- 1 file changed, 59 insertions(+), 46 deletions(-) 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)] From 9d62278bb319a121c16e649a1c4ee83c655b6b41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 11 Oct 2025 23:13:12 +0200 Subject: [PATCH 2/3] Use stack buffers for CPU exponent lowering and add MutableUInt128 --- PerfectNumbers.Core/MutableUInt128.cs | 108 +++++++++++++++++ PerfectNumbers.Core/PrimeOrderCalculator.cs | 123 +++++++++++--------- 2 files changed, 177 insertions(+), 54 deletions(-) create mode 100644 PerfectNumbers.Core/MutableUInt128.cs diff --git a/PerfectNumbers.Core/MutableUInt128.cs b/PerfectNumbers.Core/MutableUInt128.cs new file mode 100644 index 0000000..7922ec1 --- /dev/null +++ b/PerfectNumbers.Core/MutableUInt128.cs @@ -0,0 +1,108 @@ +using System; +using System.Numerics; + +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 sealed class MutableUInt128 +{ + public ulong High { get; private set; } + public ulong Low { get; private set; } + + public MutableUInt128() + { + } + + public MutableUInt128(ulong value) + { + Low = value; + } + + public MutableUInt128(ulong high, ulong low) + { + High = high; + Low = low; + } + + public MutableUInt128(UInt128 value) + { + Set(value); + } + + public void Set(UInt128 value) + { + High = (ulong)(value >> 64); + Low = (ulong)value; + } + + public void Set(ulong high, ulong low) + { + High = high; + Low = low; + } + + public UInt128 ToUInt128() => ((UInt128)High << 64) | Low; + + public void Add(ulong value) + { + ulong low = Low + value; + if (low < value) + { + High++; + } + + Low = low; + } + + 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; + } + + public void Subtract(ulong value) + { + ulong low = Low; + ulong result = low - value; + ulong borrow = result > low ? 1UL : 0UL; + + Low = result; + High -= borrow; + } + + 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; + } + + public void Multiply(ulong value) + { + ulong high = High; + ulong low = Low; + + ulong highContribution = Math.BigMul(low, value, out ulong newLow); + Math.BigMul(high, value, out ulong shiftedHigh); + ulong newHigh = shiftedHigh + highContribution; + + Low = newLow; + High = newHigh; + } + + public void MultiplyAdd(ulong multiplier, ulong addend) + { + Multiply(multiplier); + Add(addend); + } +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 1312815..7df88c5 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -314,69 +314,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) { From 60c8dbf6d6abb5f68947461981c52ac68f187efb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 11 Oct 2025 23:13:17 +0200 Subject: [PATCH 3/3] Make MutableUInt128 a struct and reuse it in Pollard Rho --- PerfectNumbers.Core/MutableUInt128.cs | 79 +++++++++++++++++---- PerfectNumbers.Core/PrimeOrderCalculator.cs | 11 ++- 2 files changed, 74 insertions(+), 16 deletions(-) diff --git a/PerfectNumbers.Core/MutableUInt128.cs b/PerfectNumbers.Core/MutableUInt128.cs index 7922ec1..1b53b11 100644 --- a/PerfectNumbers.Core/MutableUInt128.cs +++ b/PerfectNumbers.Core/MutableUInt128.cs @@ -1,5 +1,6 @@ using System; using System.Numerics; +using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; @@ -7,17 +8,14 @@ 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 sealed class MutableUInt128 +internal struct MutableUInt128 { - public ulong High { get; private set; } - public ulong Low { get; private set; } - - public MutableUInt128() - { - } + public ulong High; + public ulong Low; public MutableUInt128(ulong value) { + High = 0UL; Low = value; } @@ -29,23 +27,35 @@ internal sealed class MutableUInt128 public MutableUInt128(UInt128 value) { - Set(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; @@ -57,6 +67,7 @@ internal sealed class MutableUInt128 Low = low; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Add(in MutableUInt128 value) { ulong low = Low + value.Low; @@ -67,6 +78,7 @@ internal sealed class MutableUInt128 High = high; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Subtract(ulong value) { ulong low = Low; @@ -77,6 +89,7 @@ internal sealed class MutableUInt128 High -= borrow; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Subtract(in MutableUInt128 value) { ulong low = Low; @@ -87,22 +100,62 @@ internal sealed class MutableUInt128 High = High - value.High - borrow; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Multiply(ulong value) { - ulong high = High; - ulong low = Low; + ulong currentHigh = High; + ulong currentLow = Low; - ulong highContribution = Math.BigMul(low, value, out ulong newLow); - Math.BigMul(high, value, out ulong shiftedHigh); - ulong newHigh = shiftedHigh + highContribution; + 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 7df88c5..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 @@ -1326,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)