Merge pull request #111 from mlemanczyk/codex/integrate-primeordercalculator-into-evenperfectbitscanner

Add heuristic UInt128 prime order calculation

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

@@ -28,15 +28,18 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
     private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
     private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
     private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
+    private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
     private readonly ulong[] _smallExponents = new ulong[SampleCount];
     private readonly ulong[] _largeExponents = new ulong[SampleCount];
     private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
     private readonly ulong[] _smallCycles = new ulong[SampleCount];
     private readonly ulong[] _largeCycles = new ulong[SampleCount];
     private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
+    private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
 
     private readonly Random _random = new(113);
     private ulong? _previousPrimeOrder;
+    private UInt128? _previousWidePrimeOrder;
 
     public enum InputScale
     {
@@ -73,6 +76,10 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
             Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
 #endif
             _veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
+
+            UInt128 wideModulus = NextWideOddModulus();
+            _wideModuli[i] = wideModulus;
+            _wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
         }
     }
 
@@ -147,6 +154,21 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
         return checksum;
     }
 
+    [Benchmark]
+    public UInt128 HeuristicWideCycleComputation()
+    {
+        UInt128 checksum = UInt128.Zero;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
+            _wideCycles[i] = cycle;
+            checksum ^= cycle;
+        }
+
+        return checksum;
+    }
+
     private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
     {
         switch (Scale)
@@ -204,6 +226,17 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
         return value;
     }
 
+    private UInt128 NextWideOddModulus()
+    {
+        Span<byte> buffer = stackalloc byte[16];
+        _random.NextBytes(buffer);
+        ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+        Span<byte> highSlice = buffer.Slice(8, 8);
+        ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
+        UInt128 value = ((UInt128)high << 64) | low;
+        return value | UInt128.One;
+    }
+
     private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
     {
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -261,6 +294,33 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
         return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
     }
 
+    private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
+    {
+        if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
+        {
+            return MersenneDivisorCycles.GetCycle(modulus);
+        }
+
+#if DEBUG
+        Console.WriteLine("Trying heuristic. Wide prime order calculation");
+#endif
+        PrimeOrderCalculator.PrimeOrderResultWide orderResult = PrimeOrderCalculator.Calculate(
+            modulus,
+            _previousWidePrimeOrder,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+
+        if (orderResult.Order != UInt128.Zero)
+        {
+            _previousWidePrimeOrder = orderResult.Order;
+            return orderResult.Order;
+        }
+
+        _previousWidePrimeOrder = null;
+
+        Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
+        return MersenneDivisorCycles.GetCycle(modulus);
+    }
+
     private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
 
     private static ulong ComputeMontgomeryNPrime(ulong modulus)

PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs

@@ -1,4 +1,5 @@
 using System.Collections.Generic;
+using System.Numerics;
 using FluentAssertions;
 using Xunit;
 
@@ -142,6 +143,29 @@ public class PrimeOrderCalculatorTests
         strict.Order.Should().Be(heuristic.Order);
     }
 
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void Calculate_handles_128_bit_prime()
+    {
+        UInt128 prime = UInt128.Parse("18446744073709641691");
+
+        PrimeOrderCalculator.PrimeOrderResultWide result = PrimeOrderCalculator.Calculate(
+            prime,
+            previousOrder: null,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be((UInt128)1229782938247309446UL);
+    }
+
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void GetCycle_uses_wide_prime_heuristic()
+    {
+        UInt128 prime = UInt128.Parse("18446744073709641691");
+        MersenneDivisorCycles.GetCycle(prime).Should().Be((UInt128)1229782938247309446UL);
+    }
+
     private static ulong ComputeOrderByDoubling(ulong prime)
     {
         if (prime == 2UL)

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -126,30 +126,48 @@ public class MersenneDivisorCycles
                 return CalculateCycleLength(divisor);
 	}
 
-	public static UInt128 GetCycle(UInt128 divisor)
-	{
-		// For divisor = 2^k, cycle is 1
-		if ((divisor & (divisor - UInt128.One)) == UInt128.Zero)
-		{
-			return UInt128.One;
-		}
+    public static UInt128 GetCycle(UInt128 divisor)
+    {
+        // For divisor = 2^k, cycle is 1
+        if ((divisor & (divisor - UInt128.One)) == UInt128.Zero)
+        {
+            return UInt128.One;
+        }
 
-		// Otherwise, find order of 2 mod divisor
-                // TODO: Port this UInt128 path to the unrolled-hex cycle calculator so wide
-                // divisors stop relying on the slow shift-and-subtract loop measured in the
-                // GPU benchmarks.
-		UInt128 order = UInt128.One, pow = UInt128Numbers.Two;
-		while (pow != UInt128.One)
-		{
-			pow <<= 1;
-			if (pow >= divisor)
-				pow -= divisor;
+        if (divisor <= ulong.MaxValue)
+        {
+            ulong cycle = CalculateCycleLength((ulong)divisor);
+            return (UInt128)cycle;
+        }
 
-			order++;
-		}
+        PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
+            divisor,
+            previousOrder: null,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+        if (wideResult.Order != UInt128.Zero)
+        {
+            return wideResult.Order;
+        }
 
-		return order;
-	}
+        // Otherwise, find order of 2 mod divisor
+        // TODO: Port this UInt128 path to the unrolled-hex cycle calculator so wide
+        // divisors stop relying on the slow shift-and-subtract loop measured in the
+        // GPU benchmarks.
+        UInt128 order = UInt128.One;
+        UInt128 pow = UInt128Numbers.Two;
+        while (pow != UInt128.One)
+        {
+            pow <<= 1;
+            if (pow >= divisor)
+            {
+                pow -= divisor;
+            }
+
+            order++;
+        }
+
+        return order;
+    }
 
     public static bool TryCalculateCycleLengthForExponent(
         ulong divisor,

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -9,7 +9,7 @@ using System.Threading;
 
 namespace PerfectNumbers.Core;
 
-internal static class PrimeOrderCalculator
+internal static partial class PrimeOrderCalculator
 {
     internal enum PrimeOrderStatus
     {