Merge pull request #156 from mlemanczyk/codex/analyze-gc-pressure-in-evenperfectbitscanner

Reduce GC churn in CPU Montgomery pow2 ladder

PerfectNumbers.Core/MutableUInt128.cs

@@ -0,0 +1,161 @@
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+/// <summary>
+/// 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.
+/// </summary>
+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);
+    }
+}

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]
+    [ThreadStatic]
-	private static ArrayPool<FactorEntry>? _factorPool;
+    private static ArrayPool<FactorEntry>? _factorPool;
-	private static ArrayPool<FactorEntry> FactorPool => _factorPool ??= ArrayPool<FactorEntry>.Create();
+    private static ArrayPool<FactorEntry> FactorPool => _factorPool ??= ArrayPool<FactorEntry>.Create();
 
-	[ThreadStatic]
+    private const int FactorEntryStackThreshold = 32;
-	private static ArrayPool<ulong>? _ulongPool;
-	private static ArrayPool<ulong> ULongPool => _ulongPool ??= ArrayPool<ulong>.Create();
 
-	private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+    [ThreadStatic]
-	{
+    private static ArrayPool<ulong>? _ulongPool;
-		ArrayPool<FactorEntry> factorPool = FactorPool;
+    private static ArrayPool<ulong> ULongPool => _ulongPool ??= ArrayPool<ulong>.Create();
 
-		FactorEntry[]? tempArray = null;
+    private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
-		try
+    {
-		{
+        ArrayPool<FactorEntry> factorPool = FactorPool;
-			ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
-			int length = factors.Count;
-			tempArray = factorPool.Rent(length + 1);
-			Span<FactorEntry> buffer = tempArray;
-			factorSpan.CopyTo(buffer);
 
-			bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor);
+        FactorEntry[]? tempArray = null;
-			// bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None);
+        try
-			if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime)
+        {
-			{
+            ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
-				buffer[length] = new FactorEntry(factors.Cofactor);
+            int length = factors.Count;
-				length++;
+            int capacity = length + 1;
-			}
 
-			buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value));
+            Span<FactorEntry> 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++)
+            factorSpan.CopyTo(buffer);
-			{
-				ulong primeFactor = buffer[i].Value;
-				int exponent = buffer[i].Exponent;
-				for (int iteration = 0; iteration < exponent; iteration++)
-				{
-					if ((order % primeFactor) != 0UL)
-					{
-						break;
-					}
 
-					ulong reduced = order / primeFactor;
+            bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor);
-					if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+            // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None);
-					{
+            if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime)
-						order = reduced;
+            {
-						continue;
+                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;
+            for (int i = 0; i < buffer.Length; i++)
-		}
+            {
-		finally
+                FactorEntry entry = buffer[i];
-		{
+                ulong primeFactor = entry.Value;
-			if (tempArray is not null)
+                int exponent = entry.Exponent;
-			{
+                for (int iteration = 0; iteration < exponent; iteration++)
-				factorPool.Return(tempArray, clearArray: false);
+                {
-			}
+                    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;
+        ref MutableUInt128 polynomial = ref s_pollardRhoPolynomialBuffer;
-        return (ulong)(value % modulus);
+        polynomial.Set(x);
+        polynomial.Multiply(x);
+        polynomial.Add(c);
+        return polynomial.Mod(modulus);
     }
 
     private static ulong BinaryGcd(ulong a, ulong 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<ulong> 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<ulong>.Shared.Rent(oddPowerCount);
+            oddPowerStorage = pooledOddPowers;
+        }
+
+        Span<ulong> oddPowers = oddPowerStorage[..oddPowerCount];
+        InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers);
 
         int index = bitLength - 1;
-		while (index >= 0)
+        while (index >= 0)
-		{
+        {
-			if (((exponent >> index) & 1UL) == 0UL)
+            if (((exponent >> index) & 1UL) == 0UL)
-			{
+            {
-				result = result.MontgomeryMultiply(result, modulus, nPrime);
+                result = result.MontgomeryMultiply(result, modulus, nPrime);
-				index--;
+                index--;
-				continue;
+                continue;
-			}
+            }
 
-			int windowStart = index - windowSize + 1;
+            int windowStart = index - windowSize + 1;
-			if (windowStart < 0)
+            if (windowStart < 0)
-			{
+            {
-				windowStart = 0;
+                windowStart = 0;
-			}
+            }
 
-			while (((exponent >> windowStart) & 1UL) == 0UL)
+            while (((exponent >> windowStart) & 1UL) == 0UL)
-			{
+            {
-				windowStart++;
+                windowStart++;
-			}
+            }
 
-			int windowLength = index - windowStart + 1;
+            int windowLength = index - windowStart + 1;
-			for (int square = 0; square < windowLength; square++)
+            for (int square = 0; square < windowLength; square++)
-			{
+            {
-				result = result.MontgomeryMultiply(result, modulus, nPrime);
+                result = result.MontgomeryMultiply(result, modulus, nPrime);
-			}
+            }
 
-			ulong mask = (1UL << windowLength) - 1UL;
+            ulong mask = (1UL << windowLength) - 1UL;
-			ulong windowValue = (exponent >> windowStart) & mask;
+            ulong windowValue = (exponent >> windowStart) & mask;
-			int tableIndex = (int)((windowValue - 1UL) >> 1);
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
-			ulong multiplier = oddPowers[tableIndex];
+            ulong multiplier = oddPowers[tableIndex];
-			result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
+            result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
 
-			index = windowStart - 1;
+            index = windowStart - 1;
-		}
+        }
-		
+
-		if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold)
+        if (pooledOddPowers is not null)
-		{
+        {
-			ArrayPool<ulong>.Shared.Return(oddPowers);
+            ArrayPool<ulong>.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<ulong> destination)
     {
-        ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold
+        int oddPowerCount = destination.Length;
-            ? new ulong[oddPowerCount]
-            : ArrayPool<ulong>.Shared.Rent(oddPowerCount);
-
         if (oddPowerCount == 0)
         {
-            return oddPowers;
+            return;
         }
 
-        Span<ulong> 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)]