Merge pull request #114 from mlemanczyk/codex/implement-first-todo-from-todo.txt

Rely on Pollard Rho for GPU prime detection

PerfectNumbers.Core/DivisorCycleCache.cs

@@ -262,9 +262,39 @@ public sealed class DivisorCycleCache
         }
 
         ulong divisor = divisors[index];
+
+        if ((divisor & (divisor - 1UL)) == 0UL)
+        {
+            cycles[index] = 1UL;
+            status[index] = ByteOne;
+            pow[index] = 1UL;
+            order[index] = 1UL;
+            return;
+        }
+
         ulong currentPow = pow[index];
         ulong currentOrder = order[index];
 
+        if (divisor <= 3UL)
+        {
+            while (currentPow != 1UL)
+            {
+                currentPow += currentPow;
+                if (currentPow >= divisor)
+                {
+                    currentPow -= divisor;
+                }
+
+                currentOrder++;
+            }
+
+            cycles[index] = currentOrder;
+            status[index] = ByteOne;
+            pow[index] = currentPow;
+            order[index] = currentOrder;
+            return;
+        }
+
         do
         {
             currentPow += currentPow;

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -143,7 +143,8 @@ public class MersenneDivisorCycles
         PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
         if (wideResult.Order != UInt128.Zero)
         {
             return wideResult.Order;
@@ -192,7 +193,8 @@ public class MersenneDivisorCycles
         PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
         if (orderResult.Order != 0UL)
         {
             cycleLength = orderResult.Order;
@@ -818,15 +820,24 @@ public class MersenneDivisorCycles
 
         public static ulong CalculateCycleLength(ulong divisor)
         {
-                // For divisor = 2^k, cycle is 1
+                if (TryCalculateCycleLengthHeuristic(divisor, out ulong cycleLength))
+                        return cycleLength;
+
+                return CalculateCycleLengthFallback(divisor);
+        }
+
+        internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, out ulong cycleLength)
+        {
                 if ((divisor & (divisor - 1UL)) == 0UL)
                 {
-                        return 1UL;
+                        cycleLength = 1UL;
+                        return true;
                 }
 
                 if (divisor <= 3UL)
                 {
-                        return CalculateCycleLengthFallback(divisor);
+                        cycleLength = CalculateCycleLengthFallback(divisor);
+                        return true;
                 }
 
                 if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
@@ -834,14 +845,17 @@ public class MersenneDivisorCycles
                         PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
                                 divisor,
                                 previousOrder: null,
-                                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+                                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+                                PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
                         if (orderResult.Order != 0UL)
                         {
-                                return orderResult.Order;
+                                cycleLength = orderResult.Order;
+                                return true;
                         }
                 }
 
-                return CalculateCycleLengthFallback(divisor);
+                cycleLength = 0UL;
+                return false;
         }
 
         private static ulong CalculateCycleLengthFallback(ulong divisor)

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -23,7 +23,15 @@ internal static partial class PrimeOrderCalculator
     }
 
     public static PrimeOrderResultWide Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
+        => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
+
+    public static PrimeOrderResultWide Calculate(
+        UInt128 prime,
+        UInt128? previousOrder,
+        PrimeOrderSearchConfig config,
+        PrimeOrderHeuristicDevice device)
     {
+        using var scope = UsePow2Mode(device);
         if (prime <= ulong.MaxValue)
         {
             ulong? previous = null;
@@ -40,7 +48,7 @@ internal static partial class PrimeOrderCalculator
                 }
             }
 
-            PrimeOrderResult result = Calculate((ulong)prime, previous, config);
+            PrimeOrderResult result = Calculate((ulong)prime, previous, config, device);
             return new PrimeOrderResultWide(result.Status, result.Order);
         }
 
@@ -60,12 +68,23 @@ internal static partial class PrimeOrderCalculator
         }
 
         UInt128 phi = prime - UInt128.One;
+        using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
         PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
         if (phiFactors.Factors is null)
         {
             return FinishStrictlyWide(prime, config.Mode);
         }
 
+        return RunHeuristicPipelineWide(prime, previousOrder, config, phi, phiFactors);
+    }
+
+    private static PrimeOrderResultWide RunHeuristicPipelineWide(
+        UInt128 prime,
+        UInt128? previousOrder,
+        PrimeOrderSearchConfig config,
+        UInt128 phi,
+        PartialFactorResult128 phiFactors)
+    {
         if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, phiFactors))
         {
             return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
@@ -553,41 +572,29 @@ internal static partial class PrimeOrderCalculator
         }
 
         Dictionary<UInt128, int> counts = new(capacity: 8);
-        UInt128 remaining = value;
-        uint[] primes = PrimesGenerator.SmallPrimes;
-        ulong[] squares = PrimesGenerator.SmallPrimesPow2;
-        int primeCount = primes.Length;
         uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
+        UInt128 remaining;
 
-        for (int i = 0; i < primeCount; i++)
+        if (IsGpuPow2Allowed && value <= ulong.MaxValue)
         {
-            uint primeCandidate = primes[i];
+            Dictionary<ulong, int> narrowCounts = new(capacity: 8);
-            if (primeCandidate > limit)
+            if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining))
             {
-                break;
+                foreach (KeyValuePair<ulong, int> entry in narrowCounts)
+                {
+                    counts[(UInt128)entry.Key] = entry.Value;
                 }
 
-            UInt128 primeValue = primeCandidate;
+                remaining = (UInt128)narrowRemaining;
-            UInt128 primeSquare = (UInt128)squares[i];
-            if (primeSquare > remaining)
-            {
-                break;
             }
-
+            else
-            if ((remaining % primeValue) != UInt128.Zero)
             {
-                continue;
+                remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
             }
-
+        }
-            int exponent = 0;
+        else
-            do
         {
-                remaining /= primeValue;
+            remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
-                exponent++;
-            }
-            while ((remaining % primeValue) == UInt128.Zero);
-
-            counts[primeValue] = exponent;
         }
 
         List<UInt128> pending = new();
@@ -682,6 +689,46 @@ internal static partial class PrimeOrderCalculator
         return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index);
     }
 
+    private static UInt128 PopulateSmallPrimeFactorsCpuWide(UInt128 value, uint limit, Dictionary<UInt128, int> counts)
+    {
+        UInt128 remaining = value;
+        uint[] primes = PrimesGenerator.SmallPrimes;
+        ulong[] squares = PrimesGenerator.SmallPrimesPow2;
+        int primeCount = primes.Length;
+
+        for (int i = 0; i < primeCount; i++)
+        {
+            uint primeCandidate = primes[i];
+            if (primeCandidate > limit)
+            {
+                break;
+            }
+
+            UInt128 primeSquare = (UInt128)squares[i];
+            if (primeSquare > remaining)
+            {
+                break;
+            }
+
+            UInt128 primeValue = primeCandidate;
+            if ((remaining % primeValue) != UInt128.Zero)
+            {
+                continue;
+            }
+
+            int exponent = 0;
+            do
+            {
+                remaining /= primeValue;
+                exponent++;
+            }
+            while ((remaining % primeValue) == UInt128.Zero);
+
+            counts[primeValue] = exponent;
+        }
+
+        return remaining;
+    }
     private static bool TryPollardRhoWide(UInt128 n, Stopwatch stopwatch, long budgetTicks, out UInt128 factor)
     {
         factor = UInt128.Zero;
@@ -955,6 +1002,9 @@ internal static partial class PrimeOrderCalculator
         {
             return UInt128.Zero;
         }
+
+        if (IsGpuPow2Allowed)
+        {
             if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue)
             {
                 ulong prime64 = (ulong)modulus;
@@ -971,6 +1021,7 @@ internal static partial class PrimeOrderCalculator
             {
                 return wideRemainder;
             }
+        }
 
         return exponent.Pow2MontgomeryModWindowed(modulus);
     }

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Buffers;
 using System.Buffers.Binary;
 using System.Collections.Generic;
@@ -61,7 +62,104 @@ internal static partial class PrimeOrderCalculator
         public static PrimeOrderSearchConfig StrictDefault => new(1_000_000, 0, 0, PrimeOrderMode.Strict);
     }
 
+    internal enum PrimeOrderHeuristicDevice
+    {
+        Cpu,
+        Gpu,
+    }
+
+    [ThreadStatic]
+    private static bool s_pow2ModeInitialized;
+
+    [ThreadStatic]
+    private static bool s_allowGpuPow2;
+
+    private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
+
+
+    [ThreadStatic]
+    private static PrimeOrderHeuristicDevice s_deviceMode;
+
+    [ThreadStatic]
+    private static bool s_debugLoggingEnabled;
+
+    private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
+    private readonly struct Pow2ModeScope : IDisposable
+    {
+        private readonly bool _previousInitialized;
+        private readonly bool _previousAllow;
+        private readonly PrimeOrderHeuristicDevice _previousDevice;
+
+        public Pow2ModeScope(PrimeOrderHeuristicDevice device)
+        {
+            _previousInitialized = s_pow2ModeInitialized;
+            _previousAllow = s_allowGpuPow2;
+            _previousDevice = s_deviceMode;
+            s_pow2ModeInitialized = true;
+            s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu;
+            s_deviceMode = device;
+        }
+
+        public void Dispose()
+        {
+            s_allowGpuPow2 = _previousAllow;
+            s_pow2ModeInitialized = _previousInitialized;
+            s_deviceMode = _previousDevice;
+        }
+    }
+
+    private readonly struct DebugLoggingScope : IDisposable
+    {
+        private readonly bool _previous;
+
+        public DebugLoggingScope(bool enabled)
+        {
+            _previous = s_debugLoggingEnabled;
+            s_debugLoggingEnabled = enabled;
+        }
+
+        public void Dispose()
+        {
+            s_debugLoggingEnabled = _previous;
+        }
+    }
+
+    private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device);
+
+    private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled);
+
+    [Conditional("DEBUG")]
+    private static void DebugLog(string message)
+    {
+        if (s_debugLoggingEnabled)
+        {
+            Console.WriteLine(message);
+        }
+    }
+
+    [Conditional("DEBUG")]
+    private static void DebugLog(Func<string> messageFactory)
+    {
+        if (s_debugLoggingEnabled)
+        {
+            Console.WriteLine(messageFactory());
+        }
+    }
+
     public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config)
+        => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
+
+    public static PrimeOrderResult Calculate(
+        ulong prime,
+        ulong? previousOrder,
+        PrimeOrderSearchConfig config,
+        PrimeOrderHeuristicDevice device)
+    {
+        using var scope = UsePow2Mode(device);
+        return CalculateInternal(prime, previousOrder, config);
+    }
+
+    private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config)
     {
         if (prime <= 3UL)
         {
@@ -71,36 +169,44 @@ internal static partial class PrimeOrderCalculator
         ulong phi = prime - 1UL;
         MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(prime);
 
-#if DEBUG
+        if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult))
-        Console.WriteLine("Partial factoring φ(p)");
+        {
-#endif
+            return gpuResult;
+        }
+
+        using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
+        DebugLog("Partial factoring φ(p)");
+
         PartialFactorResult phiFactors = PartialFactor(phi, config);
 
         if (phiFactors.Factors is null)
         {
-#if DEBUG
+            DebugLog("No factors found");
-            Console.WriteLine("No factors found");
-#endif
             return FinishStrictly(prime, config.Mode);
         }
 
-#if DEBUG
+        return RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors);
-        Console.WriteLine("Trying special max check");
+    }
-#endif
+
+    private static PrimeOrderResult RunHeuristicPipeline(
+        ulong prime,
+        ulong? previousOrder,
+        PrimeOrderSearchConfig config,
+        in MontgomeryDivisorData divisorData,
+        ulong phi,
+        PartialFactorResult phiFactors)
+    {
+        DebugLog("Trying special max check");
         if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData))
         {
             return new PrimeOrderResult(PrimeOrderStatus.Found, phi);
         }
 
-#if DEBUG
+        DebugLog("Initializing starting order");
-        Console.WriteLine("Initializing starting order");
-#endif
         ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData);
         candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData);
 
-#if DEBUG
+        DebugLog("Trying to confirm order");
-        Console.WriteLine("Trying to confirm order");
-#endif
         if (TryConfirmOrder(prime, candidateOrder, divisorData, config))
         {
             return new PrimeOrderResult(PrimeOrderStatus.Found, candidateOrder);
@@ -116,9 +222,7 @@ internal static partial class PrimeOrderCalculator
             return new PrimeOrderResult(PrimeOrderStatus.Found, order);
         }
 
-#if DEBUG
+        DebugLog("Heuristic unresolved, finishing strictly");
-        Console.WriteLine("Heuristic unresolved, finishing strictly");
-#endif
         return FinishStrictly(prime, config.Mode);
     }
 
@@ -222,17 +326,13 @@ internal static partial class PrimeOrderCalculator
         }
 
         // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`.
-#if DEBUG
+        DebugLog("Verifying a^order ≡ 1 (mod p)");
-        Console.WriteLine("Verifying a^order ≡ 1 (mod p)");
-#endif
         if (!Pow2EqualsOne(order, prime, divisorData))
         {
             return false;
         }
 
-#if DEBUG
+        DebugLog("Partial factoring order");
-        Console.WriteLine("Partial factoring order");
-#endif
 
         // TODO: Do we do partial factoring of order multiple times?
         PartialFactorResult factorization = PartialFactor(order, config);
@@ -245,16 +345,12 @@ internal static partial class PrimeOrderCalculator
         {
             if (factorization.Cofactor <= 1UL)
             {
-#if DEBUG
+            DebugLog("Cofactor <= 1. No factors found");
-                Console.WriteLine("Cofactor <= 1. No factors found");
-#endif
                 return false;
             }
 
             // TODO: Use Open.Numerics.Primality for this final check once it's available.
-#if DEBUG
+            DebugLog("Cofactor > 1. Testing primality of cofactor");
-            Console.WriteLine("Cofactor > 1. Testing primality of cofactor");
-#endif
             bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
             // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
             if (!isPrime)
@@ -262,16 +358,12 @@ internal static partial class PrimeOrderCalculator
                 return false;
             }
 
-#if DEBUG
+            DebugLog("Adding cofactor as prime factor");
-            Console.WriteLine("Adding cofactor as prime factor");
-#endif
             factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
         }
 
         ReadOnlySpan<FactorEntry> span = factorization.Factors;
-#if DEBUG
+        DebugLog("Verifying prime-power reductions");
-        Console.WriteLine("Verifying prime-power reductions");
-#endif
         int length = factorization.Count;
         for (int i = 0; i < length; i++)
         {
@@ -311,9 +403,7 @@ internal static partial class PrimeOrderCalculator
         }
 
         // TODO: Do we do partial factoring of order multiple times?
-#if DEBUG
+        DebugLog("Trying heuristic. Partial factoring order");
-        Console.WriteLine("Trying heuristic. Partial factoring order");
-#endif
         PartialFactorResult orderFactors = PartialFactor(order, config);
         if (orderFactors.Factors is null)
         {
@@ -340,18 +430,14 @@ internal static partial class PrimeOrderCalculator
         int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
         List<ulong> candidates = new(capacity);
         FactorEntry[] factorArray = orderFactors.Factors!;
-#if DEBUG
+        DebugLog("Building candidates list");
-        Console.WriteLine("Building candidates list");
-#endif
         BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
         if (candidates.Count == 0)
         {
             return false;
         }
 
-#if DEBUG
+        DebugLog("Sorting candidates");
-        Console.WriteLine("Sorting candidates");
-#endif
         SortCandidates(prime, previousOrder, candidates);
 
         int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
@@ -360,9 +446,7 @@ internal static partial class PrimeOrderCalculator
         bool allowGpuBatch = true;
         Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
 
-#if DEBUG
+        DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
-        Console.WriteLine($"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
-#endif
         int index = 0;
         const int MaxGpuBatchSize = 256;
         const int StackGpuBatchSize = 32;
@@ -384,7 +468,7 @@ internal static partial class PrimeOrderCalculator
             bool gpuStackRemainders = false;
             GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
 
-            if (allowGpuBatch)
+            if (allowGpuBatch && IsGpuPow2Allowed)
             {
                 if (batchSize <= StackGpuBatchSize)
                 {
@@ -461,9 +545,7 @@ internal static partial class PrimeOrderCalculator
             index += batchSize;
         }
 
-#if DEBUG
+        DebugLog("No candidate confirmed");
-        Console.WriteLine("No candidate confirmed");
-#endif
         return false;
     }
 
@@ -690,12 +772,15 @@ internal static partial class PrimeOrderCalculator
     }
 
     private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
+    {
+        if (IsGpuPow2Allowed)
         {
             GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder);
             if (status == GpuPow2ModStatus.Success)
             {
                 return remainder == 1UL;
             }
+        }
 
         return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
     }
@@ -707,20 +792,261 @@ internal static partial class PrimeOrderCalculator
             return PartialFactorResult.Empty;
         }
 
-        Dictionary<ulong, int> counts = new(capacity: 8);
+        const int FactorSlotCount = GpuSmallPrimeFactorSlots;
+        Span<ulong> primeSlots = stackalloc ulong[FactorSlotCount];
+        Span<int> exponentSlots = stackalloc int[FactorSlotCount];
+        primeSlots.Clear();
+        exponentSlots.Clear();
+
+        int factorCount = 0;
+        Dictionary<ulong, int>? counts = null;
+        bool useDictionary = false;
+
+        uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
+        ulong remaining;
+
+        bool gpuFactored = false;
+        if (IsGpuHeuristicDevice)
+        {
+            gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
+        }
+        else
+        {
+            remaining = value;
+        }
+
+        if (!gpuFactored)
+        {
+            primeSlots.Clear();
+            exponentSlots.Clear();
+            if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
+            {
+                useDictionary = true;
+                counts = new Dictionary<ulong, int>(capacity: 8);
+                remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
+            }
+        }
+
+        List<ulong> pending = new();
+        if (remaining > 1UL)
+        {
+            pending.Add(remaining);
+        }
+
+        if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
+        {
+            DebugLog("Processing pending composites with Pollard's Rho");
+            Stopwatch stopwatch = Stopwatch.StartNew();
+            long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
+            Stack<ulong> stack = new();
+            stack.Push(remaining);
+            pending.Clear();
+
+            while (stack.Count > 0)
+            {
+                ulong composite = stack.Pop();
+                if (composite <= 1UL)
+                {
+                    continue;
+                }
+
+                bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
+                if (isPrime)
+                {
+                    AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
+                    continue;
+                }
+
+                if (stopwatch.ElapsedTicks > budgetTicks)
+                {
+                    pending.Add(composite);
+                    continue;
+                }
+
+                if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor))
+                {
+                    pending.Add(composite);
+                    continue;
+                }
+
+                ulong quotient = composite / factor;
+                stack.Push(factor);
+                stack.Push(quotient);
+            }
+        }
+
+        ulong cofactor = 1UL;
+        int pendingCount = pending.Count;
+        for (int i = 0; i < pendingCount; i++)
+        {
+            ulong composite = pending[i];
+            bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
+            if (isPrime)
+            {
+                AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
+            }
+            else if (composite > 1UL)
+            {
+                cofactor = checked(cofactor * composite);
+            }
+        }
+
+        if (useDictionary)
+        {
+            if ((counts is null || counts.Count == 0) && cofactor == value)
+            {
+                DebugLog("cofactor is the same as value, no factors found");
+                return new PartialFactorResult(null, value, false, 0);
+            }
+
+            if (counts is null || counts.Count == 0)
+            {
+                return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
+            }
+
+            FactorEntry[] factors = ArrayPool<FactorEntry>.Shared.Rent(counts.Count);
+            int index = 0;
+            DebugLog(() => $"Collecting {counts.Count} prime factors");
+            foreach (KeyValuePair<ulong, int> entry in counts)
+            {
+                factors[index] = new FactorEntry(entry.Key, entry.Value);
+                index++;
+            }
+
+            Span<FactorEntry> span = factors.AsSpan(0, index);
+            span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+
+            FactorEntry[] resultArray = new FactorEntry[index];
+            span.CopyTo(resultArray);
+            ArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
+
+            bool fullyFactored = cofactor == 1UL;
+            return new PartialFactorResult(resultArray, cofactor, fullyFactored, index);
+        }
+
+        int actualCount = 0;
+        for (int i = 0; i < factorCount; i++)
+        {
+            if (primeSlots[i] != 0UL && exponentSlots[i] != 0)
+            {
+                actualCount++;
+            }
+        }
+
+        if (actualCount == 0)
+        {
+            if (cofactor == value)
+            {
+                DebugLog("cofactor is the same as value, no factors found");
+                return new PartialFactorResult(null, value, false, 0);
+            }
+
+            return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
+        }
+
+        FactorEntry[] array = new FactorEntry[actualCount];
+        int arrayIndex = 0;
+        DebugLog(() => $"Collecting {actualCount} prime factors");
+        for (int i = 0; i < factorCount; i++)
+        {
+            ulong primeValue = primeSlots[i];
+            int exponentValue = exponentSlots[i];
+            if (primeValue == 0UL || exponentValue == 0)
+            {
+                continue;
+            }
+
+            array[arrayIndex] = new FactorEntry(primeValue, exponentValue);
+            arrayIndex++;
+        }
+
+        Array.Sort(array, static (a, b) => a.Value.CompareTo(b.Value));
+        bool fullyFactoredArray = cofactor == 1UL;
+        return new PartialFactorResult(array, cofactor, fullyFactoredArray, array.Length);
+    }
+
+    private const int GpuSmallPrimeFactorSlots = 64;
+
+    private static bool TryPopulateSmallPrimeFactorsCpu(
+        ulong value,
+        uint limit,
+        Span<ulong> primeTargets,
+        Span<int> exponentTargets,
+        out int factorCount,
+        out ulong remaining)
+    {
+        factorCount = 0;
+        remaining = value;
+
+        if (value <= 1UL)
+        {
+            return true;
+        }
+
+        int capacity = Math.Min(primeTargets.Length, exponentTargets.Length);
+        if (capacity == 0)
+        {
+            return false;
+        }
+
+        uint[] primes = PrimesGenerator.SmallPrimes;
+        ulong[] squares = PrimesGenerator.SmallPrimesPow2;
+        int primeCount = primes.Length;
+        ulong remainingLocal = value;
+        uint effectiveLimit = limit == 0 ? uint.MaxValue : limit;
+
+        for (int i = 0; i < primeCount && remainingLocal > 1UL; i++)
+        {
+            uint primeCandidate = primes[i];
+            if (primeCandidate > effectiveLimit)
+            {
+                break;
+            }
+
+            ulong primeSquare = squares[i];
+            if (primeSquare != 0UL && primeSquare > remainingLocal)
+            {
+                break;
+            }
+
+            ulong primeValue = primeCandidate;
+            if ((remainingLocal % primeValue) != 0UL)
+            {
+                continue;
+            }
+
+            int exponent = 0;
+            do
+            {
+                remainingLocal /= primeValue;
+                exponent++;
+            }
+            while ((remainingLocal % primeValue) == 0UL);
+
+            if (factorCount >= capacity)
+            {
+                return false;
+            }
+
+            primeTargets[factorCount] = primeValue;
+            exponentTargets[factorCount] = exponent;
+            factorCount++;
+        }
+
+        remaining = remainingLocal;
+        return true;
+    }
+
+    private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary<ulong, int> counts)
+    {
         ulong remaining = value;
         uint[] primes = PrimesGenerator.SmallPrimes;
         ulong[] squares = PrimesGenerator.SmallPrimesPow2;
         int primeCount = primes.Length;
-        uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
 
         for (int i = 0; i < primeCount; i++)
         {
             uint primeCandidate = primes[i];
-
-#if DEBUG
-            Console.WriteLine($"Trying small prime {primeCandidate}");
-#endif
             if (primeCandidate > limit)
             {
                 break;
@@ -748,107 +1074,37 @@ internal static partial class PrimeOrderCalculator
             counts[primeValue] = exponent;
         }
 
-        List<ulong> pending = new();
+        return remaining;
-        if (remaining > 1UL)
-        {
-            pending.Add(remaining);
     }
 
-        if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
+    private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out ulong remaining)
     {
-#if DEBUG
+        Span<ulong> primeBuffer = stackalloc ulong[GpuSmallPrimeFactorSlots];
-            Console.WriteLine("Processing pending composites with Pollard's Rho");
+        Span<int> exponentBuffer = stackalloc int[GpuSmallPrimeFactorSlots];
-#endif
+        primeBuffer.Clear();
-            Stopwatch stopwatch = Stopwatch.StartNew();
+        exponentBuffer.Clear();
-            long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
-            Stack<ulong> stack = new();
-            stack.Push(remaining);
-            pending.Clear();
 
-            while (stack.Count > 0)
+        if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _))
         {
-                ulong composite = stack.Pop();
+            remaining = value;
-                if (composite == 1UL)
+            return false;
+        }
+
+        remaining = gpuRemaining;
+        for (int i = 0; i < factorCount; i++)
+        {
+            ulong primeValue = primeBuffer[i];
+            int exponent = exponentBuffer[i];
+            if (primeValue == 0UL || exponent == 0)
             {
                 continue;
             }
 
-                bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
+            counts[primeValue] = exponent;
-                // if (PrimeTester.IsPrimeInternal(composite, CancellationToken.None))
-                if (isPrime)
-                {
-                    AddFactor(counts, composite, 1);
-                    continue;
         }
 
-                if (stopwatch.ElapsedTicks > budgetTicks)
+        return true;
-                {
-                    pending.Add(composite);
-                    continue;
     }
-
-                if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor))
-                {
-                    pending.Add(composite);
-                    continue;
-                }
-
-                ulong quotient = composite / factor;
-                stack.Push(factor);
-                stack.Push(quotient);
-            }
-        }
-
-        ulong cofactor = 1UL;
-        int pendingCount = pending.Count;
-#if DEBUG
-        Console.WriteLine($"Processing {pendingCount} pending composites with Open.Numeric.Primes");
-#endif
-        for (int i = 0; i < pendingCount; i++)
-        {
-            ulong composite = pending[i];
-            bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
-            // if (PrimeTester.IsPrimeInternal(composite, CancellationToken.None))
-            if (isPrime)
-            {
-                AddFactor(counts, composite, 1);
-            }
-            else
-            {
-                cofactor = checked(cofactor * composite);
-            }
-        }
-
-        if (counts.Count == 0 && cofactor == value)
-        {
-#if DEBUG
-            Console.WriteLine("cofactor is the same as value, no factors found");
-#endif
-            return new PartialFactorResult(null, value, false, 0);
-        }
-
-        FactorEntry[] factors = ArrayPool<FactorEntry>.Shared.Rent(counts.Count);
-        int index = 0;
-#if DEBUG
-        Console.WriteLine($"Collecting {counts.Count} prime factors");
-#endif
-        foreach (KeyValuePair<ulong, int> entry in counts)
-        {
-            factors[index] = new FactorEntry(entry.Key, entry.Value);
-            index++;
-        }
-
-        Span<FactorEntry> span = factors.AsSpan(0, index);
-        span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
-
-        FactorEntry[] resultArray = new FactorEntry[index];
-        span.CopyTo(resultArray);
-        ArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
-
-        bool fullyFactored = cofactor == 1UL;
-        return new PartialFactorResult(resultArray, cofactor, fullyFactored, index);
-    }
-
     private static bool TryPollardRho(ulong n, Stopwatch stopwatch, long budgetTicks, out ulong factor)
     {
         factor = 0UL;
@@ -946,6 +1202,84 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
+    private static bool TryAppendFactor(Span<ulong> primes, Span<int> exponents, ref int count, ulong prime, int exponent)
+    {
+        if (exponent <= 0)
+        {
+            return true;
+        }
+
+        for (int i = 0; i < count; i++)
+        {
+            if (primes[i] == prime)
+            {
+                exponents[i] += exponent;
+                return true;
+            }
+        }
+
+        if (count >= primes.Length)
+        {
+            return false;
+        }
+
+        primes[count] = prime;
+        exponents[count] = exponent;
+        count++;
+        return true;
+    }
+
+    private static void AddFactorToCollector(
+        ref bool useDictionary,
+        ref Dictionary<ulong, int>? counts,
+        Span<ulong> primes,
+        Span<int> exponents,
+        ref int count,
+        ulong prime,
+        int exponent)
+    {
+        if (exponent <= 0)
+        {
+            return;
+        }
+
+        if (useDictionary)
+        {
+            AddFactor(counts!, prime, exponent);
+            return;
+        }
+
+        if (TryAppendFactor(primes, exponents, ref count, prime, exponent))
+        {
+            return;
+        }
+
+        counts ??= new Dictionary<ulong, int>(Math.Max(count, 8));
+        CopyFactorsToDictionary(primes, exponents, count, counts);
+        count = 0;
+        useDictionary = true;
+        AddFactor(counts, prime, exponent);
+    }
+
+    private static void CopyFactorsToDictionary(
+        ReadOnlySpan<ulong> primes,
+        ReadOnlySpan<int> exponents,
+        int count,
+        Dictionary<ulong, int> target)
+    {
+        for (int i = 0; i < count; i++)
+        {
+            ulong prime = primes[i];
+            int exponent = exponents[i];
+            if (prime == 0UL || exponent == 0)
+            {
+                continue;
+            }
+
+            target[prime] = exponent;
+        }
+    }
+
     private static ulong CalculateByDoubling(ulong prime)
     {
         ulong order = 1UL;

todo.txt

@@ -1,2 +1 @@
-- Mirror the complete CPU heuristic cycle algorithm on the GPU path so both implementations stay aligned.
 - Extend the GPU numeric type probe to cover proposed custom types such as GpuUInt256 and GpuUInt512 once their definitions are available so we can plan the wide-order heuristics on the device.

tools/GpuNumericTypeProbe.csx

@@ -1,5 +1,6 @@
 #r "nuget: ILGPU, 1.5.3"
 #r "nuget: PeterO.Numbers, 1.8.1"
+#r "nuget: Open.Numeric.Primes, 4.0.4"
 
 using System;
 using System.Numerics;
@@ -7,6 +8,7 @@ using ILGPU;
 using ILGPU.Runtime;
 using ILGPU.Runtime.CPU;
 using PeterO.Numbers;
+using Open.Numeric.Primes;
 
 var context = Context.CreateDefault();
 var accelerator = context.CreateCPUAccelerator(0);
@@ -20,7 +22,8 @@ try
     [
         ProbeBigInteger(accelerator),
         ProbeEInteger(accelerator),
-        ProbeERational(accelerator)
+        ProbeERational(accelerator),
+        ProbeOpenNumericPrimesIsPrime(accelerator)
     ];
 
     foreach (ProbeResult result in results)
@@ -112,6 +115,31 @@ static ProbeResult ProbeERational(Accelerator accelerator)
     }
 }
 
+static ProbeResult ProbeOpenNumericPrimesIsPrime(Accelerator accelerator)
+{
+    try
+    {
+        var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<int, Stride1D.Dense>>(OpenNumericPrimesIsPrimeKernel);
+        using var buffer = accelerator.Allocate1D<int>(4);
+        kernel(new Index1D((int)buffer.Length), buffer.View);
+        accelerator.Synchronize();
+        return ProbeResult.Success(
+            "Open.Numeric.Primes.Prime.Numbers.IsPrime",
+            "IsPrime compiles and runs on the CPU accelerator kernel.");
+    }
+    catch (Exception ex)
+    {
+        return ProbeResult.Failure("Open.Numeric.Primes.Prime.Numbers.IsPrime", ex);
+    }
+}
+
+static void OpenNumericPrimesIsPrimeKernel(Index1D index, ArrayView1D<int, Stride1D.Dense> data)
+{
+    ulong value = (ulong)(index + 2);
+    bool isPrime = Prime.Numbers.IsPrime(value);
+    data[index] = isPrime ? 1 : 0;
+}
+
 static void BigIntegerKernel(Index1D index, ArrayView1D<int, Stride1D.Dense> data)
 {
     BigInteger value = new BigInteger(data[index]);