From 2bf9ff233c91875ff72fc21a70591938377c9a86 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Sun, 21 Dec 2025 18:02:59 +0100 Subject: [PATCH] PrimeOrderCalculator converted to a template. --- .../Pow2ModBenchmarks.cs | 24 +- .../Gpu/PrimeOrderGpuHeuristicsTests.cs | 2 +- .../PrimeOrderCalculatorTests.cs | 32 +- .../Gpu/PrimeOrderConstants.cs | 4 + .../Gpu/PrimeOrderGpuHeuristics.cs | 2 +- .../MersenneDivisorCyclesTemplate.cs | 61 +- PerfectNumbers.Core/PerfectNumberConstants.cs | 1 + .../PerfectNumbers.Core.csproj | 1 + PerfectNumbers.Core/PrimeNumberConstants.cs | 6 + .../PrimeOrderCalculator.Config.cs | 39 - .../PrimeOrderCalculator.Cpu.cs | 1909 ---------- .../PrimeOrderCalculator.Gpu.cs | 1041 ------ .../PrimeOrderCalculator.Hybrid.cs | 651 ---- .../PrimeOrderCalculator.Shared.cs | 631 ---- .../PrimeOrderCalculatorConfig.cs | 34 + .../PrimeOrderCalculatorTemplate.cs | 3273 +++++++++++++++++ PerfectNumbers.Core/ThreadStaticPools.cs | 25 + PerfectNumbers.Core/UInt128Numbers.cs | 1 + PerfectNumbers.Core/UnboundedTaskScheduler.cs | 2 +- 19 files changed, 3407 insertions(+), 4332 deletions(-) create mode 100644 PerfectNumbers.Core/PrimeNumberConstants.cs delete mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Config.cs delete mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs delete mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs delete mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs delete mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs create mode 100644 PerfectNumbers.Core/PrimeOrderCalculatorConfig.cs create mode 100644 PerfectNumbers.Core/PrimeOrderCalculatorTemplate.cs diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs index 0149d29..7c3d39b 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs @@ -619,11 +619,11 @@ public class Pow2ModBenchmarks return MersenneDivisorCyclesCpu.CalculateCycleLength(modulus, divisorData); } - ulong order = PrimeOrderCalculator.CalculateCpu( + ulong order = PrimeOrderCalculatorCpu.Calculate( modulus, previousPrimeOrder, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (order != 0UL) { @@ -642,10 +642,10 @@ public class Pow2ModBenchmarks return MersenneDivisorCyclesCpu.GetCycle(modulus); } - UInt128 order = PrimeOrderCalculator.CalculateCpu( + UInt128 order = PrimeOrderCalculatorCpu.Calculate( modulus, previousWidePrimeOrder, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (order != UInt128.Zero) { @@ -689,11 +689,11 @@ public class Pow2ModBenchmarks #if DEBUG Console.WriteLine("Trying heuristic. Prime order calculation"); #endif - ulong order = PrimeOrderCalculator.CalculateCpu( + ulong order = PrimeOrderCalculatorCpu.Calculate( modulus, _previousPrimeOrder, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (order != 0UL) { @@ -718,12 +718,12 @@ public class Pow2ModBenchmarks #if DEBUG Console.WriteLine("Trying heuristic. Prime order calculation"); #endif - ulong order = PrimeOrderCalculator.CalculateGpu( + ulong order = PrimeOrderCalculatorGpu.Calculate( gpu, modulus, _previousPrimeOrder, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (order != 0UL) { @@ -747,10 +747,10 @@ public class Pow2ModBenchmarks #if DEBUG Console.WriteLine("Trying heuristic. Wide prime order calculation"); #endif - UInt128 order = PrimeOrderCalculator.CalculateCpu( + UInt128 order = PrimeOrderCalculatorCpu.Calculate( modulus, _previousWidePrimeOrder, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (order != UInt128.Zero) { @@ -774,11 +774,11 @@ public class Pow2ModBenchmarks #if DEBUG Console.WriteLine("Trying heuristic. Wide prime order calculation"); #endif - UInt128 order = PrimeOrderCalculator.CalculateGpu( + UInt128 order = PrimeOrderCalculatorGpu.Calculate( gpu, modulus, _previousWidePrimeOrder, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (order != UInt128.Zero) { diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs index a9ebc5d..7381817 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs @@ -398,7 +398,7 @@ public class PrimeOrderGpuHeuristicsTests PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); - ulong result = PrimeOrderCalculator.CalculateGpu(gpu, prime, null, divisorData, PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + ulong result = PrimeOrderCalculatorGpu.Calculate(gpu, prime, null, divisorData, PrimeOrderCalculatorConfig.HeuristicDefault); result.Should().Be(3UL); } diff --git a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs index 9d0f1bd..2330d14 100644 --- a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs @@ -45,12 +45,12 @@ public class PrimeOrderCalculatorTests ulong expected = ComputeOrderByDoubling(prime); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - ulong result = PrimeOrderCalculator.CalculateGpu( + ulong result = PrimeOrderCalculatorGpu.Calculate( gpu, prime, previousOrder: null, MontgomeryDivisorData.FromModulus(prime), - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorAccelerator.Return(gpu); @@ -65,19 +65,19 @@ public class PrimeOrderCalculatorTests ulong previousOrder = ComputeOrderByDoubling(previousPrime); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - ulong heuristic = PrimeOrderCalculator.CalculateGpu( + ulong heuristic = PrimeOrderCalculatorGpu.Calculate( gpu, prime, previousOrder, MontgomeryDivisorData.FromModulus(prime), - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); - ulong strict = PrimeOrderCalculator.CalculateGpu( + ulong strict = PrimeOrderCalculatorGpu.Calculate( gpu, prime, previousOrder: null, MontgomeryDivisorData.FromModulus(prime), - PrimeOrderCalculator.PrimeOrderCalculatorConfig.StrictDefault); + PrimeOrderCalculatorConfig.StrictDefault); PrimeOrderCalculatorAccelerator.Return(gpu); @@ -91,12 +91,12 @@ public class PrimeOrderCalculatorTests public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder) { var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - ulong result = PrimeOrderCalculator.CalculateGpu( + ulong result = PrimeOrderCalculatorGpu.Calculate( gpu, prime, previousOrder: null, MontgomeryDivisorData.FromModulus(prime), - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorAccelerator.Return(gpu); result.Should().Be(expectedOrder); @@ -106,14 +106,14 @@ public class PrimeOrderCalculatorTests [Trait("Category", "Fast")] public void Calculate_heuristic_mode_returns_unresolved_when_phi_cannot_be_factored() { - var config = new PrimeOrderCalculator.PrimeOrderCalculatorConfig( + var config = new PrimeOrderCalculatorConfig( smallFactorLimit: 1, pollardRhoMilliseconds: 0, maxPowChecks: 8, strictMode: false); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - ulong result = PrimeOrderCalculator.CalculateGpu( + ulong result = PrimeOrderCalculatorGpu.Calculate( gpu, prime: 13UL, previousOrder: null, @@ -128,7 +128,7 @@ public class PrimeOrderCalculatorTests [Trait("Category", "Fast")] public void Calculate_heuristic_mode_attempts_candidates_before_falling_back_to_strict_mode() { - var heuristicConfig = new PrimeOrderCalculator.PrimeOrderCalculatorConfig( + var heuristicConfig = new PrimeOrderCalculatorConfig( smallFactorLimit: 2, pollardRhoMilliseconds: 0, maxPowChecks: 1, @@ -138,7 +138,7 @@ public class PrimeOrderCalculatorTests try { MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL); - ulong heuristic = PrimeOrderCalculator.CalculateGpu( + ulong heuristic = PrimeOrderCalculatorGpu.Calculate( gpu, prime: 239UL, previousOrder: null, @@ -147,13 +147,13 @@ public class PrimeOrderCalculatorTests heuristic.Should().Be(ComputeOrderByDoubling(239UL)); - var strictConfig = new PrimeOrderCalculator.PrimeOrderCalculatorConfig( + var strictConfig = new PrimeOrderCalculatorConfig( smallFactorLimit: 2, pollardRhoMilliseconds: 0, maxPowChecks: 1, strictMode: true); - ulong strict = PrimeOrderCalculator.CalculateGpu( + ulong strict = PrimeOrderCalculatorGpu.Calculate( gpu, prime: 239UL, previousOrder: null, @@ -175,11 +175,11 @@ public class PrimeOrderCalculatorTests UInt128 prime = UInt128.Parse("18446744073709641691"); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - UInt128 result = PrimeOrderCalculator.CalculateGpu( + UInt128 result = PrimeOrderCalculatorGpu.Calculate( gpu, prime, previousOrder: null, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorAccelerator.Return(gpu); result.Should().Be((UInt128)1229782938247309446UL); diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderConstants.cs b/PerfectNumbers.Core/Gpu/PrimeOrderConstants.cs index 1bb8678..487e7de 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderConstants.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderConstants.cs @@ -10,5 +10,9 @@ public static class PrimeOrderConstants public const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); public const int Pow2WindowSizeBits = 8; public const int WideStackThreshold = 12; + public const int MaxCompositeSlots = 256; + public const int ExponentSlotLimit = 256; + public const int PrimeSlotLimit = 256; + public const int ScratchSlots = 64; } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index e580cbc..5294afc 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -253,7 +253,7 @@ internal static partial class PrimeOrderGpuHeuristics PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong? previousOrder, - PrimeOrderCalculator.PrimeOrderCalculatorConfig config, + PrimeOrderCalculatorConfig config, in MontgomeryDivisorData divisorData, out ulong order) { diff --git a/PerfectNumbers.Core/MersenneDivisorCyclesTemplate.cs b/PerfectNumbers.Core/MersenneDivisorCyclesTemplate.cs index 18198b9..419a90f 100644 --- a/PerfectNumbers.Core/MersenneDivisorCyclesTemplate.cs +++ b/PerfectNumbers.Core/MersenneDivisorCyclesTemplate.cs @@ -259,22 +259,23 @@ public struct MersenneDivisorCyclesTemplate } #if DEVICE_GPU - UInt128 wideOrder = PrimeOrderCalculator.CalculateGpu( + UInt128 wideOrder = PrimeOrderCalculatorGpu.Calculate( gpu, divisor, previousOrder: null, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); #elif DEVICE_HYBRID - UInt128 wideOrder = PrimeOrderCalculator.CalculateHybrid( + UInt128 wideOrder = PrimeOrderCalculatorHybrid.Calculate( gpu, divisor, previousOrder: null, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + MontgomeryDivisorData.FromModulus((ulong)divisor), + PrimeOrderCalculatorConfig.HeuristicDefault); #else - UInt128 wideOrder = PrimeOrderCalculator.CalculateCpu( + UInt128 wideOrder = PrimeOrderCalculatorCpu.Calculate( divisor, previousOrder: null, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); #endif if (wideOrder != UInt128.Zero) { @@ -329,25 +330,25 @@ public struct MersenneDivisorCyclesTemplate // } #if DEVICE_GPU - ulong computedOrder = PrimeOrderCalculator.CalculateGpu( + ulong computedOrder = PrimeOrderCalculatorGpu.Calculate( gpu, divisor, previousOrder: null, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); #elif DEVICE_HYBRID - ulong computedOrder = PrimeOrderCalculator.CalculateHybrid( + ulong computedOrder = PrimeOrderCalculatorHybrid.Calculate( gpu, divisor, previousOrder: null, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); #else - ulong computedOrder = PrimeOrderCalculator.CalculateCpu( + ulong computedOrder = PrimeOrderCalculatorCpu.Calculate( divisor, previousOrder: null, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); #endif if (computedOrder != 0UL) { @@ -1011,14 +1012,14 @@ public struct MersenneDivisorCyclesTemplate } #if DEVICE_GPU - if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeGpu(gpu, divisor)) - { - ulong computedOrder = PrimeOrderCalculator.CalculateGpu( - gpu, - divisor, - previousOrder: null, - divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeGpu(gpu, divisor)) + { + ulong computedOrder = PrimeOrderCalculatorGpu.Calculate( + gpu, + divisor, + previousOrder: null, + divisorData, + PrimeOrderCalculatorConfig.HeuristicDefault); if (computedOrder != 0UL) { cycleLength = computedOrder; @@ -1026,14 +1027,14 @@ public struct MersenneDivisorCyclesTemplate } } #elif DEVICE_HYBRID - if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor)) - { - ulong computedOrder = PrimeOrderCalculator.CalculateHybrid( - gpu, - divisor, - previousOrder: null, - divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor)) + { + ulong computedOrder = PrimeOrderCalculatorHybrid.Calculate( + gpu, + divisor, + previousOrder: null, + divisorData, + PrimeOrderCalculatorConfig.HeuristicDefault); if (computedOrder != 0UL) { cycleLength = computedOrder; @@ -1043,11 +1044,11 @@ public struct MersenneDivisorCyclesTemplate #else if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor)) { - ulong computedOrder = PrimeOrderCalculator.CalculateCpu( + ulong computedOrder = PrimeOrderCalculatorCpu.Calculate( divisor, previousOrder: null, divisorData, - PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); + PrimeOrderCalculatorConfig.HeuristicDefault); if (computedOrder != 0UL) { cycleLength = computedOrder; diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 45f619d..6731b8d 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -23,6 +23,7 @@ public static class PerfectNumberConstants public const int ThreadsByAccelerator = 1; public const int ByDivisorStateSaveInterval = 200_000; public const string ByDivisorStateDirectory = "Checks"; + public const int PollardRhoMilliseconds = 100; // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal // divisor-cycle datasets we generate offline. // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune diff --git a/PerfectNumbers.Core/PerfectNumbers.Core.csproj b/PerfectNumbers.Core/PerfectNumbers.Core.csproj index f4a97b8..7f393c0 100644 --- a/PerfectNumbers.Core/PerfectNumbers.Core.csproj +++ b/PerfectNumbers.Core/PerfectNumbers.Core.csproj @@ -5,6 +5,7 @@ enable enable 12.0 + $(DefineConstants);DEVICE_CPU true $(IntermediateOutputPath)generated diff --git a/PerfectNumbers.Core/PrimeNumberConstants.cs b/PerfectNumbers.Core/PrimeNumberConstants.cs new file mode 100644 index 0000000..7df204c --- /dev/null +++ b/PerfectNumbers.Core/PrimeNumberConstants.cs @@ -0,0 +1,6 @@ +namespace PerfectNumbers.Core; + +public static class PrimeNumberConstants +{ + public const ulong GpuSmallPrimeBound = 1_000_000UL; +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Config.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Config.cs deleted file mode 100644 index 9b73069..0000000 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Config.cs +++ /dev/null @@ -1,39 +0,0 @@ -using System.Runtime.CompilerServices; - -namespace PerfectNumbers.Core; - -internal static partial class PrimeOrderCalculator -{ - internal readonly struct PrimeOrderCalculatorConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, bool strictMode) - { - private static PrimeOrderCalculatorConfig _heuristicDefault = CreateHeuristicDefault(UnboundedTaskScheduler.ConfiguredThreadCount); - - public static PrimeOrderCalculatorConfig HeuristicDefault - { - [MethodImpl(MethodImplOptions.AggressiveInlining)] - get => _heuristicDefault; - } - - public static readonly PrimeOrderCalculatorConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, true); - - public static void ConfigureHeuristicDefault(int threadCount) - { - _heuristicDefault = CreateHeuristicDefault(threadCount); - } - - public static PrimeOrderCalculatorConfig CreateHeuristicDefault(int threadCount) - { - int normalizedThreadCount = threadCount < 1 ? 1 : threadCount; - return new(smallFactorLimit: 100_000, pollardRhoMilliseconds: normalizedThreadCount * 24, maxPowChecks: 24, false); - } - - public readonly uint SmallFactorLimit = smallFactorLimit; - public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; - public readonly int MaxPowChecks = maxPowChecks; - public readonly int MaxPowChecksCapacity = maxPowChecks << 2; - public readonly bool StrictMode = strictMode; - } - - // TODO: Remove branching on the CPU / GPU paths - // TODO: Split big non-static functions into smaller / extract static code to limit JIT time -} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs deleted file mode 100644 index f0673f7..0000000 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ /dev/null @@ -1,1909 +0,0 @@ -using System.Buffers; -using System.Collections.Concurrent; -using System.Diagnostics; -using System.Numerics; -using System.Runtime.CompilerServices; -using System.Runtime.InteropServices; -using PerfectNumbers.Core.Cpu; -using PerfectNumbers.Core.Gpu; - -namespace PerfectNumbers.Core; - -internal static partial class PrimeOrderCalculator -{ - public static ulong CalculateCpu( - ulong prime, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config) - { - // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. - if (prime <= 3UL) - { - return prime == 3UL ? 2UL : 1UL; - } - - ulong phi = prime - 1UL; - - PartialFactorResult phiFactors = PartialFactorCpu(phi, divisorData, config); - - ulong result; - if (phiFactors.Factors is null) - { - result = CalculateByFactorizationCpu(prime, divisorData); - - phiFactors.Dispose(); - return result; - } - - result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); - phiFactors.Dispose(); - return result; - } - - public static UInt128 CalculateCpu( - in UInt128 prime, - in UInt128? previousOrder, - in PrimeOrderCalculatorConfig config) - { - UInt128 result; - if (prime <= ulong.MaxValue) - { - ulong? previous = null; - if (previousOrder.HasValue) - { - UInt128 previousValue = previousOrder.Value; - if (previousValue <= ulong.MaxValue) - { - previous = (ulong)previousValue; - } - else - { - previous = ulong.MaxValue; - } - } - - ulong prime64 = (ulong)prime; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime64); - ulong order64 = CalculateCpu(prime64, previous, divisorData, config); - result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; - } - else - { - result = CalculateWideInternalCpu(prime, previousOrder, MontgomeryDivisorData.Empty, config); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 CalculateWideInternalCpu(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) - { - if (prime <= UInt128.One) - { - return UInt128.One; - } - - if (prime == (UInt128)3UL) - { - return (UInt128)2UL; - } - - UInt128 phi = prime - UInt128.One; - PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); - UInt128 result; - if (phiFactors.Factors is null) - { - result = FinishStrictlyWideCpu(prime, divisorData); - } - else - { - result = RunHeuristicPipelineWideCpu(prime, previousOrder, divisorData, config, phi, phiFactors); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong RunHeuristicPipelineCpu( - ulong prime, - ulong? previousOrder, - in PrimeOrderCalculatorConfig config, - in MontgomeryDivisorData divisorData, - ulong phi, - PartialFactorResult phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) - { - return phi; - } - - ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); - candidateOrder = ExponentLoweringCpu(candidateOrder, phiFactors, divisorData); - - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) - { - return candidateOrder; - } - - if (config.StrictMode) - { - orderFactors?.Dispose(); - return CalculateByFactorizationCpu(prime, divisorData); - } - - if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, orderFactors, out ulong order)) - { - return order; - } - - return candidateOrder; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 RunHeuristicPipelineWideCpu( - in UInt128 prime, - in UInt128? previousOrder, - in MontgomeryDivisorData? divisorData, - in PrimeOrderCalculatorConfig config, - in UInt128 phi, - in PartialFactorResult128 phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxWideCpu(phi, prime, divisorData, phiFactors)) - { - return phi; - } - - UInt128 candidateOrder = InitializeStartingOrderWideCpu(prime, phi, divisorData); - candidateOrder = ExponentLoweringWideCpu(candidateOrder, prime, divisorData, phiFactors); - - if (TryConfirmOrderWideCpu(prime, candidateOrder, divisorData, config)) - { - return candidateOrder; - } - - if (config.StrictMode) - { - return FinishStrictlyWideCpu(prime, divisorData); - } - - if (TryHeuristicFinishWideCpu(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) - { - return order; - } - - return candidateOrder; - } - - private static readonly ConcurrentDictionary _dictionaryStats = new(20_480, 64); - private static void PrintStats(int count) - { - _dictionaryStats.AddOrUpdate(count, 1UL, (index, item) => item + 1UL); - Console.WriteLine("Current stats:"); - foreach(var (key, value) in _dictionaryStats) - { - Console.WriteLine($"Count {key} = {value}"); - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 FinishStrictlyWideCpu(in UInt128 prime, in MontgomeryDivisorData? divisorData) - { - UInt128 phi = prime - UInt128.One; - FixedCapacityStack> uInt128IntDictionaryPool = ThreadStaticPools.UInt128IntDictionaryPool; - Dictionary counts = uInt128IntDictionaryPool.Rent(capacity: 8); - FactorCompletelyWide(phi, counts); - int entryCount = counts.Count; - PrintStats(entryCount); - if (entryCount == 0) - { - uInt128IntDictionaryPool.Return(counts); - return phi; - } - - FixedCapacityStack>> keyValuePairUInt128IntegerPool = ThreadStaticPools.KeyValuePairUInt128IntegerPool; - var entries = keyValuePairUInt128IntegerPool.Rent(entryCount); - entries.AddRange(counts); - uInt128IntDictionaryPool.Return(counts); - - entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); - - UInt128 order = phi; - for (int i = 0; i < entryCount; i++) - { - UInt128 primeFactor = entries[i].Key; - int exponent = entries[i].Value; - for (int iteration = 0; iteration < exponent; iteration++) - { - if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero) - { - UInt128 candidate = order / primeFactor; - if (Pow2ModWideCpu(candidate, prime, divisorData) == UInt128.One) - { - order = candidate; - } - else - { - break; - } - } - } - } - - keyValuePairUInt128IntegerPool.Return(entries); - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxWideCpu(in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - for (int i = 0; i < length; i++) - { - UInt128 factor = factorSpan[i].Value; - UInt128 reduced = phi / factor; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 InitializeStartingOrderWideCpu(in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData? divisorData) - { - UInt128 order = phi; - UInt128 mod8 = prime & (UInt128)7UL; - if (mod8 == UInt128.One || mod8 == (UInt128)7UL) - { - UInt128 half = phi >> 1; - if (Pow2ModWideCpu(half, prime, divisorData) == UInt128.One) - { - order = half; - } - } - - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 ExponentLoweringWideCpu(UInt128 order, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) - { - ArrayPool pool = ThreadStaticPools.FactorEntry128Pool; - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - // TODO: This should never trigger from production code - check - // if (length == 0) - // { - // return order; - // } - FactorEntry128[] tempArray = pool.Rent(length + 1); - Span buffer = tempArray.AsSpan(0, length); - - factorSpan.CopyTo(buffer); - if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor)) - { - buffer[length] = new FactorEntry128(factors.Cofactor, 1); - length++; - } - - buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value)); - - for (int i = 0; i < length; i++) - { - UInt128 primeFactor = buffer[i].Value; - int exponent = buffer[i].Exponent; - for (int iteration = 0; iteration < exponent; iteration++) - { - if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero) - { - UInt128 reduced = order / primeFactor; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) - { - order = reduced; - continue; - } - } - - break; - } - } - - pool.Return(tempArray, clearArray: false); - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmOrderWideCpu(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) - { - if (order == UInt128.Zero) - { - return false; - } - - if (Pow2ModWideCpu(order, prime, divisorData) != UInt128.One) - { - return false; - } - - PartialFactorResult128 factorization = PartialFactorWide(order, config); - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(factorization.Cofactor)) - { - return false; - } - - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - UInt128 primeFactor = span[i].Value; - UInt128 reduced = order; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) - { - break; - } - - reduced /= primeFactor; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - } - - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryHeuristicFinishWideCpu( - in UInt128 prime, - in UInt128 order, - in UInt128? previousOrder, - in MontgomeryDivisorData? divisorData, - in PrimeOrderCalculatorConfig config, - out UInt128 result) - { - result = UInt128.Zero; - if (order <= UInt128.One) - { - return false; - } - - PartialFactorResult128 orderFactors = PartialFactorWide(order, config); - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(orderFactors.Cofactor)) - { - return false; - } - - orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - FixedCapacityStack> uInt128ListPool = ThreadStaticPools.UInt128ListPool; - List candidates = uInt128ListPool.Rent(capacity); - FactorEntry128[] factorArray = orderFactors.Factors!; - - BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); - - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - uInt128ListPool.Return(candidates); - return false; - } - - SortCandidatesWide(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - - for (int i = 0; i < candidateCount; i++) - { - if (powUsed >= powBudget && powBudget > 0) - { - break; - } - - UInt128 candidate = candidates[i]; - powUsed++; - - if (Pow2ModWideCpu(candidate, prime, divisorData) != UInt128.One) - { - continue; - } - - if (!TryConfirmCandidateWideCpu(prime, candidate, config, ref powUsed, powBudget, divisorData)) - { - continue; - } - - result = candidate; - uInt128ListPool.Return(candidates); - return true; - } - - uInt128ListPool.Return(candidates); - return false; - } - - // private static ulong specialMaxHits; - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - int length = factors.Count; - // The phi factorization on the scanner path always yields at least one entry for phi >= 2, - // preventing the zero-length case from occurring outside synthetic tests. - // if (length == 0) - // { - // return true; - // } - - ReadOnlySpan factorSpan = new(factors.Factors, 0, length); - Span stackBuffer = stackalloc ulong[length]; - - return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool EvaluateSpecialMaxCandidatesCpu(Span buffer, ReadOnlySpan factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData) - { - int factorCount = factors.Length; - for (int i = 0; i < factorCount; i++) - { - ulong factor = factors[i]; - - // The partial factor pipeline never feeds zero or oversized factors while scanning candidate orders. - // if (factor == 0UL || factor > phi) - // { - // continue; - // } - - // if (reduced == 0UL) - // { - // continue; - // } - - buffer[i] = phi / factor; - } - - buffer.Sort(); - - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - - if (stepper.InitializeCpuIsUnity(buffer[0])) - { - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - return false; - } - - for (int i = 1; i < factorCount; i++) - { - if (stepper.ComputeNextIsUnity(buffer[i])) - { - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - return false; - } - } - - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) - { - ulong order = phi; - if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) - { - ulong half = phi >> 1; - if (Pow2EqualsOneCpu(half, divisorData)) - { - order = half; - } - } - - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static ulong ExponentLoweringCpu(ulong order, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - int factorCount = factors.Count; - int newFactorCount = factorCount + 1; - Span factorsSpan = new(factors.Factors, 0, newFactorCount); - Span exponentsSpan = new(factors.Exponents, 0, newFactorCount); - - if (!factors.FullyFactored && factors.Cofactor > 1UL && factors.CofactorIsPrime) - { - factorsSpan[factorCount] = factors.Cofactor; - exponentsSpan[factorCount] = 1; - factors.Count = factorCount = newFactorCount; - factors.Cofactor = 1UL; - factors.FullyFactored = true; - factors.CofactorIsPrime = false; - } - - SortFactorExponentSpans(factorsSpan, exponentsSpan, factorCount); - - ExponentRemainderStepperCpu stepper = factors.ExponentRemainderStepper; - - Span stackCandidates = new(factors.StackCandidates); - Span stackEvaluations = new(factors.StackEvaluations); - - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factorsSpan[i]; - int exponent = exponentsSpan[i]; - // Factor exponents produced by Pollard-Rho and trial division are always positive in the scanner flow. - // if (exponent <= 0) - // { - // continue; - // } - - ProcessExponentLoweringPrime(stackCandidates[..exponent], stackEvaluations[..exponent], ref order, primeFactor, exponent, ref stepper); - } - - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static void SortFactorExponentSpans(Span factorsSpan, Span exponentsSpan, int length) - { - if (length <= 1) - { - return; - } - - for (int i = 1; i < length; i++) - { - int current = i; - int previous = current - 1; - while (current > 0 && factorsSpan[current] < factorsSpan[previous]) - { - (factorsSpan[current], factorsSpan[previous]) = (factorsSpan[previous], factorsSpan[current]); - (exponentsSpan[current], exponentsSpan[previous]) = (exponentsSpan[previous], exponentsSpan[current]); - current--; - previous--; - } - } - } - - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void ProcessExponentLoweringPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepperCpu stepper) - { - ulong working = order; - int actual = 0; - while (actual < exponent) - { - if (working < primeFactor) - { - break; - } - - ulong reduced = working / primeFactor; - if (reduced == 0UL) - { - break; - } - - candidateBuffer[actual] = reduced; - working = reduced; - actual++; - } - - if (actual == 0) - { - return; - } - - Span candidates = candidateBuffer[..actual]; - Span evaluations = evaluationBuffer[..actual]; - - stepper.Reset(); - - int last = actual - 1; - evaluations[last] = stepper.InitializeCpuIsUnity(candidates[last]); - for (int j = last - 1; j >= 0; j--) - { - evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); - } - - for (int j = 0; j < actual; j++) - { - if (!evaluations[j]) - { - break; - } - - order = candidates[j]; - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmOrderCpu( - ulong prime, - ulong order, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - out PartialFactorResult? reusableFactorization) - { - reusableFactorization = null; - - ArgumentOutOfRangeException.ThrowIfEqual(order, 0UL); - if (!Pow2EqualsOneCpu(order, divisorData)) - { - return false; - } - - PartialFactorResult factorization = PartialFactorCpu(order, divisorData, config); - if (factorization.Factors is null) - { - factorization.Dispose(); - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) - { - factorization.Dispose(); - return false; - } - - factorization.WithAdditionalPrime(factorization.Cofactor); - } - - if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) - { - reusableFactorization = factorization; - return false; - } - - factorization.Dispose(); - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool ValidateOrderAgainstFactors( - ulong prime, - ulong order, - in MontgomeryDivisorData divisorData, - PartialFactorResult factorization) - { - ReadOnlySpan factorsSpan = factorization.Factors!; - ReadOnlySpan exponentsSpan = factorization.Exponents!; - - int length = factorization.Count; - - const int StackExponentCapacity = 32; - - Span stackBuffer = stackalloc ulong[StackExponentCapacity]; - - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool violates = false; - - for (int i = 0; i < length; i++) - { - ulong primeFactor = factorsSpan[i]; - int exponent = exponentsSpan[i]; - - // Factor exponents emitted by TryConfirmCandidateCpu stay positive for production workloads. - // if (exponent <= 0) - // { - // continue; - // } - - if (ValidateOrderForFactor(stackBuffer[..exponent], primeFactor, exponent, order, ref stepper)) - { - violates = true; - break; - } - } - - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - - return !violates; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool ValidateOrderForFactor(Span buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepperCpu stepper) - { - ulong working = order; - int actual = 0; - while (actual < exponent) - { - if (working < primeFactor) - { - break; - } - - ulong reduced = working / primeFactor; - if (reduced == 0UL) - { - break; - } - - buffer[actual] = reduced; - working = reduced; - actual++; - } - - if (actual == 0) - { - return false; - } - - Span candidates = buffer[..actual]; - - stepper.Reset(); - - int last = actual - 1; - if (stepper.InitializeCpuIsUnity(candidates[last])) - { - return true; - } - - for (int j = last - 1; j >= 0; j--) - { - if (stepper.ComputeNextIsUnity(candidates[j])) - { - return true; - } - } - - return false; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryHeuristicFinishCpu( - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - PartialFactorResult? cachedOrderFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - cachedOrderFactors?.Dispose(); - return false; - } - - // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorCpu(order, divisorData, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL) - { - return false; - } - - if (!orderFactors.CofactorIsPrime) - { - return false; - } - - orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); - candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - ThreadStaticPools.ReturnUlongList(candidates); - return false; - } - - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - - ulong candidate = candidateSpan[0]; - powUsed++; - - bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); - if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - int index = 1; - while (index < candidateCount && powUsed < powBudget) - { - int budgetRemaining = powBudget - powUsed; - if (budgetRemaining <= 0) - { - break; - } - - candidate = candidateSpan[index++]; - powUsed++; - - equalsOne = powStepper.ComputeNextIsUnity(candidate); - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - - return false; - } - finally - { - orderFactors.Dispose(); - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) - { - ulong previous = previousOrder ?? 0UL; - bool hasPrevious = previousOrder.HasValue; - ulong threshold1 = prime >> 3; - ulong threshold2 = prime >> 2; - ulong threshold3 = (prime * 3UL) >> 3; - int previousGroup = hasPrevious ? GetGroup(previous, threshold1, threshold2, threshold3) : 1; - - candidates.Sort((x, y) => CompareCandidates(x, y, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3)); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static int CompareCandidates(ulong left, ulong right, ulong previous, int previousGroup, bool hasPrevious, ulong threshold1, ulong threshold2, ulong threshold3) - { - int leftGroup = GetGroup(left, threshold1, threshold2, threshold3); - int rightGroup = GetGroup(right, threshold1, threshold2, threshold3); - if (leftGroup == 0) - { - return rightGroup == 0 ? 0 : 1; - } - - if (rightGroup == 0) - { - return -1; - } - - bool leftIsGe = !hasPrevious || left >= previous; - bool rightIsGe = !hasPrevious || right >= previous; - int rightPrimary = ComputePrimary(rightGroup, rightIsGe, previousGroup); - - int compare = ComputePrimary(leftGroup, leftIsGe, previousGroup); - compare = compare.CompareTo(rightPrimary); - if (compare != 0) - { - return compare; - } - - long leftSecondary; - long leftTertiary; - ulong tempReference; - if (leftGroup == 3) - { - leftSecondary = leftTertiary = -(long)left; - } - else - { - tempReference = hasPrevious ? (left > previous ? left - previous : previous - left) : left; - leftSecondary = (long)tempReference; - leftTertiary = (long)left; - } - - long rightSecondary; - long rightTertiary; - if (rightGroup == 3) - { - rightSecondary = rightTertiary = -(long)right; - } - else - { - tempReference = hasPrevious ? (right > previous ? right - previous : previous - right) : right; - rightSecondary = (long)tempReference; - rightTertiary = (long)right; - } - - compare = leftSecondary.CompareTo(rightSecondary); - if (compare != 0) - { - return compare; - } - - return leftTertiary.CompareTo(rightTertiary); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static int ComputePrimary(int group, bool isGe, int previousGroup) - { - int groupOffset = group switch - { - 1 => 0, - 2 => 2, - 3 => 4, - _ => 6, - }; - - if (group == previousGroup) - { - if (group == 3) - { - return groupOffset + (isGe ? 0 : 3); - } - - return groupOffset + (isGe ? 0 : 1); - } - - return groupOffset + (isGe ? 0 : 1); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static int GetGroup(ulong value, ulong threshold1, ulong threshold2, ulong threshold3) - { - if (value <= threshold1) - { - return 1; - } - - if (value <= threshold2) - { - return 2; - } - - if (value <= threshold3) - { - return 3; - } - - return 0; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static void BuildCandidates(ulong order, ulong[] factors, int[] exponents, int count, List candidates, int limit) - { - if (count == 0) - { - return; - } - - Array.Sort(factors, exponents, 0, count); - BuildCandidatesRecursive(order, factors, exponents, 0, 1UL, candidates, limit); - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void BuildCandidatesRecursive(ulong order, in ulong[] factors, in int[] exponents, int index, ulong divisorProduct, List candidates, int limit) - { - if (candidates.Count >= limit) - { - return; - } - - if (index >= factors.Length) - { - if (divisorProduct == 1UL || divisorProduct == order) - { - return; - } - - ulong candidate = order / divisorProduct; - if (candidate > 1UL && candidate < order) - { - candidates.Add(candidate); - } - - return; - } - - // FactorEntry factor = factors[index]; - // int factorExponent = factor.Exponent; - // ulong primeFactor = factor.Value; - - ulong factor = factors[index]; - int factorExponent = exponents[index]; - - ulong contribution = 1UL; - ulong contributionLimit = factorExponent == 0 ? 0UL : order / factor; - for (int exponent = 0; exponent <= factorExponent; exponent++) - { - ulong nextDivisor = divisorProduct * contribution; - if (nextDivisor > order) - { - break; - } - - BuildCandidatesRecursive(order, factors, exponents, index + 1, nextDivisor, candidates, limit); - if (candidates.Count >= limit) - { - return; - } - - if (exponent == factorExponent) - { - break; - } - - if (contribution > contributionLimit) - { - break; - } - - contribution *= factor; - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) - { - PartialFactorResult factorization = PartialFactorCpu(candidate, divisorData, config); - - if (factorization.Factors is null) - { - factorization.Dispose(); - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) - { - factorization.Dispose(); - return false; - } - - factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan factorsSpan = factorization.Factors; - ReadOnlySpan exponentsSpan = factorization.Exponents; - int length = factorization.Count; - - const int StackExponentCapacity = 32; - const int ExponentHardLimit = 256; - - Span stackBuffer = stackalloc ulong[StackExponentCapacity]; - ulong[]? heapCandidateArray = null; - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool violates = false; - - for (int i = 0; i < length; i++) - { - ulong primeFactor = factorsSpan[i]; - int exponent = exponentsSpan[i]; - - // Factorization entries for candidate confirmation are always positive here. - // if (exponent <= 0) - // { - // continue; - // } - - if (exponent > ExponentHardLimit) - { - throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); - } - - if (exponent <= StackExponentCapacity) - { - if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) - { - violates = true; - break; - } - - continue; - } - - heapCandidateArray ??= FixedCapacityPools.ExclusiveUlongArray.Rent(ExponentHardLimit); - if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) - { - violates = true; - break; - } - } - - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - - if (heapCandidateArray is not null) - { - FixedCapacityPools.ExclusiveUlongArray.Return(heapCandidateArray); - } - - factorization.Dispose(); - return !violates; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool CheckCandidateViolation(Span buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepperCpu stepper) - { - ulong working = candidate; - int actual = 0; - while (actual < exponent) - { - if (working < primeFactor) - { - break; - } - - ulong reduced = working / primeFactor; - if (reduced == 0UL) - { - break; - } - - buffer[actual] = reduced; - working = reduced; - actual++; - } - - if (actual == 0) - { - return false; - } - - Span candidates = buffer[..actual]; - - stepper.Reset(); - int last = actual - 1; - - bool enforceBudget = powBudget > 0; - if (!enforceBudget) - { - powUsed++; - if (stepper.InitializeCpuIsUnity(candidates[last])) - { - return true; - } - - for (int j = last - 1; j >= 0; j--) - { - powUsed++; - if (stepper.ComputeNextIsUnity(candidates[j])) - { - return true; - } - } - - return false; - } - - if (powUsed >= powBudget) - { - return true; - } - - powUsed++; - if (stepper.InitializeCpuIsUnity(candidates[last])) - { - return true; - } - - for (int j = last - 1; j >= 0; j--) - { - if (powUsed >= powBudget) - { - return true; - } - - powUsed++; - if (stepper.ComputeNextIsUnity(candidates[j])) - { - return true; - } - } - - return false; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static bool Pow2EqualsOneCpu(ulong exponent, in MontgomeryDivisorData divisorData) - => exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorCpu(ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - // stackalloc is faster than pooling - - var result = PartialFactorResult.Rent(divisorData); - ulong[] factors = result.Factors; - Span primeSlots = new(factors); - int[] exponents = result.Exponents; - Span exponentSlots = new(exponents); - - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); - - List pending = result.PendingFactors; - pending.Clear(); - - PopulateSmallPrimeFactorsCpu( - value, - config.SmallFactorLimit, - primeSlots, - exponentSlots, - out int factorCount, - out value); - - bool cofactorContainsComposite; - if (value > 1UL) - { - if (config.PollardRhoMilliseconds > 0) - { - FixedCapacityStack compositeStack = result.CompositeStack; - compositeStack.Clear(); - compositeStack.Push(value); - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - - CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out cofactorContainsComposite); - - if (cofactorContainsComposite) - { - while (compositeStack.Count > 0) - { - pending.Add(new(compositeStack.Pop(), knownComposite: false)); - } - } - } - else - { - pending.Add(new (value, knownComposite: false)); - } - } - - ulong cofactor = 1UL; - bool isPrime; - int pendingCount = pending.Count; - cofactorContainsComposite = false; - - for (int index = 0; index < pendingCount; index++) - { - PartialFactorPendingEntry entry = pending[index]; - // Reuse "value" parameter for testing if a number is composite to limit registry pressure. - value = entry.Value; - - if (entry.KnownComposite) - { - cofactor = checked(cofactor * value); - cofactorContainsComposite = true; - continue; - } - - if (!entry.HasKnownPrimality) - { - isPrime = PrimeTesterByLastDigit.IsPrimeCpu(value); - - entry.WithPrimality(isPrime); - pending[index] = entry; - } - - if (entry.IsPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, value); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * value); - cofactorContainsComposite = true; - } - } - - if (cofactorContainsComposite || cofactor <= 1UL) - { - isPrime = false; - } - else - { - isPrime = PrimeTesterByLastDigit.IsPrimeCpu(cofactor); - } - - // This will never happen in production code. We'll always get at least 1 factor - // if (factorCount == 0) - // { - // (...) - // } - - Array.Sort(factors, exponents, 0, factorCount); - - // Check if it was it factored completely - - result.Cofactor = cofactor; - result.FullyFactored = cofactor == 1UL; - result.Count = factorCount; - result.CofactorIsPrime = isPrime; - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CollectFactorsCpu(Span primeSlots, Span exponentSlots, ref int factorCount, List pending, FixedCapacityStack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) - { - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(composite); - - if (isPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); - continue; - } - - pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; - if (pollardRhoDeadlineReached) - { - pending.Add(new (composite, knownComposite: true)); - continue; - } - - if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(new (composite, knownComposite: true)); - continue; - } - - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - - pollardRhoDeadlineReached = false; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void PopulateSmallPrimeFactorsCpu( - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining) - { - factorCount = 0; - if (value <= 1UL) - { - remaining = value; - return; - } - - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; - - for (int i = 0; i < primeCount && value > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > limit || squares[i] > value) - { - break; - } - - // primeSquare will never == 0 in production code - // if (primeSquare != 0 && primeSquare > remainingLocal) - ulong primeValue = primeCandidate; - // primeCandidate will never equal 0 in production code - // if (primeCandidate == 0UL) - // { - // continue; - // } - - if (value != 0UL && value < primeValue) - { - continue; - } - - // This is cheap early exclusion of the candidate. We skip the candidate here in the most cases. - ulong quotient = value / primeValue; - if (value - (quotient * primeValue) != 0UL) - { - continue; - } - - int exponent = ExtractSmallPrimeExponent(ref value, primeValue, quotient); - - primeTargets[factorCount] = primeValue; - exponentTargets[factorCount] = exponent; - factorCount++; - } - - remaining = value; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - // Unroll the initial divisions because small exponents dominate in practice. Mostly we get exponent = 0-6 from this method. - private static int ExtractSmallPrimeExponent(ref ulong value, ulong primeValue, ulong quotient) - { - ulong dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - value = dividend; - return 1; - } - - dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - value = dividend; - return 2; - } - - dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - value = dividend; - return 3; - } - - dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - value = dividend; - return 4; - } - - dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - value = dividend; - return 5; - } - - dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - value = dividend; - return 6; - } - - int exponent = 7; - while (true) - { - dividend = quotient; - quotient = dividend / primeValue; - if (dividend - (quotient * primeValue) != 0UL) - { - break; - } - - exponent++; - } - - value = dividend; - return exponent; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static long CreateDeadlineTimestamp(int milliseconds) - { - if (milliseconds <= 0) - { - return long.MaxValue; - } - - long startTimestamp = Stopwatch.GetTimestamp(); - long deadline = ConvertMillisecondsToStopwatchTicks(milliseconds); - deadline += startTimestamp; - if (deadline < startTimestamp) - { - return long.MaxValue; - } - - return deadline; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static long ConvertMillisecondsToStopwatchTicks(int milliseconds) - { - long frequency = Stopwatch.Frequency; - long baseTicks = (frequency / 1000L) * milliseconds; - long remainder = (frequency % 1000L) * milliseconds; - return baseTicks + (remainder / 1000L); - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryPollardRhoCpu(ulong n, long deadlineTimestamp, out ulong factor) - { - factor = 0UL; - if ((n & 1UL) == 0UL) - { - factor = 2UL; - return true; - } - - long timestamp; // reused for deadline checks. - while (true) - { - ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } - - x = AdvancePolynomialCpu(x, c, n); - y = AdvancePolynomialTwiceCpu(y, c, n); - ulong diff = x > y ? x - y : y - x; - d = diff.BinaryGcd(n); - } - - if (d != n) - { - factor = d; - return true; - } - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static ulong AdvancePolynomialCpu(ulong x, ulong c, ulong modulus) - { - UInt128 value = (UInt128)x * x + c; - return (ulong)value.ReduceCycleRemainder(modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static ulong AdvancePolynomialTwiceCpu(ulong x, ulong c, ulong modulus) - { - UInt128 value = (UInt128)x * x + c; - x = (ulong)value.ReduceCycleRemainder(modulus); - value = (UInt128)x * x + c; - return (ulong)value.ReduceCycleRemainder(modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static ulong BinaryGcd(ulong a, ulong b) - { - if (a == 0UL) - { - return b; - } - - if (b == 0UL) - { - return a; - } - - int shift = BitOperations.TrailingZeroCount(a | b); - a >>= BitOperations.TrailingZeroCount(a); - - while (true) - { - b >>= BitOperations.TrailingZeroCount(b); - if (a > b) - { - (a, b) = (b, a); - } - - b -= a; - if (b == 0UL) - { - return a << shift; - } - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static void AddFactor(Dictionary counts, ulong prime, int exponent) - { - if (counts.TryGetValue(prime, out int existing)) - { - counts[prime] = existing + exponent; - } - else - { - counts[prime] = exponent; - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void AddFactorToCollector( - Span primes, - Span exponents, - ref int count, - ulong prime) - { - for (int i = 0; i < count; i++) - { - if (primes[i] == prime) - { - exponents[i]++; - return; - } - } - - primes[count] = prime; - exponents[count] = 1; - count++; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) - { - ulong phi = prime - 1UL; - Dictionary counts = ThreadStaticPools.RentFactorCountDictionary(); - FactorCompletelyCpu(phi, counts, false); - PrintStats(counts.Count); - if (counts.Count == 0) - { - ThreadStaticPools.ReturnFactorCountDictionary(counts); - return phi; - } - - KeyValuePair[] entries = [.. counts]; - Array.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); - - ulong order = phi; - int entryCount = entries.Length; - - for (int i = 0; i < entryCount; i++) - { - ulong primeFactor = entries[i].Key; - int exponent = entries[i].Value; - for (int iteration = 0; iteration < exponent; iteration++) - { - ulong remainder = order.ReduceCycleRemainder(primeFactor); - if (remainder != 0UL) - { - break; - } - - // TODO: Benchmark performance with the division vs reusing remainder - ulong candidate = order - primeFactor * remainder; - // ulong candidate = order / primeFactor; - - if (candidate.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL) - { - order = candidate; - continue; - } - - break; - } - } - - ThreadStaticPools.ReturnFactorCountDictionary(counts); - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void FactorCompletelyCpu(ulong value, Dictionary counts, bool knownComposite) - { - if (value <= 1UL) - { - return; - } - - bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(value); - if (!knownComposite && isPrime) - { - AddFactor(counts, value, 1); - return; - } - - ulong factor = PollardRhoStrict(value); - ulong quotient = value / factor; - - isPrime = PrimeTesterByLastDigit.IsPrimeCpu(factor); - if (isPrime) - { - int exponent = 1; - ulong remaining = quotient; - while (remaining.ReduceCycleRemainder(factor) == 0UL) - { - remaining /= factor; - exponent++; - } - - AddFactor(counts, factor, exponent); - - if (remaining > 1UL) - { - FactorCompletelyCpu(remaining, counts, knownComposite: false); - } - - return; - } - - FactorCompletelyCpu(factor, counts, knownComposite: true); - - if (quotient == factor) - { - isPrime = false; - } - else - { - isPrime = PrimeTesterByLastDigit.IsPrimeCpu(quotient); - } - - if (isPrime) - { - AddFactor(counts, quotient, 1); - } - else - { - FactorCompletelyCpu(quotient, counts, knownComposite: true); - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong PollardRhoStrict(ulong n) - { - if ((n & 1UL) == 0UL) - { - return 2UL; - } - - while (true) - { - ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - x = AdvancePolynomialCpu(x, c, n); - y = AdvancePolynomialTwiceCpu(y, c, n); - ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, n); - } - - if (d != n) - { - return d; - } - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 Pow2ModWideCpu(in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData? divisorData) - { - if (modulus == UInt128.One) - { - return UInt128.Zero; - } - - return exponent.Pow2MontgomeryModWindowed(modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmCandidateWideCpu(in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData? divisorData) - { - PartialFactorResult128 factorization = PartialFactorWide(candidate, config); - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(factorization.Cofactor)) - { - return false; - } - - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - UInt128 primeFactor = span[i].Value; - UInt128 reduced = candidate; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) - { - break; - } - - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) - { - return false; - } - - powUsed++; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - } - - return true; - } - -} - - diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs deleted file mode 100644 index 8380216..0000000 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ /dev/null @@ -1,1041 +0,0 @@ -using System.Buffers; -using System.Diagnostics; -using System.Runtime.CompilerServices; -using System.Runtime.InteropServices; -using ILGPU; -using ILGPU.Runtime; -using PerfectNumbers.Core.Gpu; -using PerfectNumbers.Core.Gpu.Accelerators; - -namespace PerfectNumbers.Core; - -internal static partial class PrimeOrderCalculator -{ - private const int GpuSmallPrimeFactorSlots = 64; - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong CalculateByFactorizationGpu(PrimeOrderCalculatorAccelerator gpu, ulong prime, in MontgomeryDivisorData divisorData) - { - ulong phi = prime - 1UL; - Dictionary counts = gpu.Pow2ModEntriesToTestOnHost; - counts.Clear(); - - FactorCompletelyCpu(phi, counts, false); - if (counts.Count == 0) - { - return phi; - } - - int entryCount = counts.Count; - Span> entries = stackalloc KeyValuePair[entryCount]; - - int stored = 0; - foreach (var entry in counts) - { - entries[stored++] = entry; - } - - entries.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); - - ulong order = phi; - var acceleratorIndex = gpu.AcceleratorIndex; - var kernel = gpu.CheckFactorsKernel; - - gpu.EnsureCapacity(entryCount, 1); - var pow2ModEntriesToTestOnDeviceView = gpu.Pow2ModEntriesToTestOnDeviceView; - - AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); - pow2ModEntriesToTestOnDeviceView.CopyFromCPU(stream, entries); - - kernel.Launch(stream, 1, entryCount, phi, pow2ModEntriesToTestOnDeviceView, divisorData.Modulus, divisorData.NPrime, divisorData.MontgomeryOne, divisorData.MontgomeryTwo, divisorData.MontgomeryTwoSquared, gpu.OutputUlongView); - - gpu.OutputUlongView.CopyToCPU(stream, ref order, 1); - stream.Synchronize(); - - AcceleratorStreamPool.Return(acceleratorIndex, stream); - return order; - } - - public static ulong CalculateGpu( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config) - { - // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. - if (prime <= 3UL) - { - return prime == 3UL ? 2UL : 1UL; - } - - ulong phi = prime - 1UL; - - if (PrimeOrderGpuHeuristics.TryCalculateOrder(gpu, prime, previousOrder, config, divisorData, out ulong gpuOrder)) - { - return gpuOrder; - } - - PartialFactorResult phiFactors = PartialFactorGpu(gpu, phi, divisorData, config); - - ulong result; - if (phiFactors.Factors is null) - { - result = CalculateByFactorizationCpu(prime, divisorData); - - phiFactors.Dispose(); - return result; - } - - result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); - phiFactors.Dispose(); - return result; - } - - public static UInt128 CalculateGpu( - PrimeOrderCalculatorAccelerator gpu, - in UInt128 prime, - in UInt128? previousOrder, - in PrimeOrderCalculatorConfig config) - { - MontgomeryDivisorData divisorData; - UInt128 result; - if (prime <= ulong.MaxValue) - { - ulong? previous = null; - if (previousOrder.HasValue) - { - UInt128 previousValue = previousOrder.Value; - if (previousValue <= ulong.MaxValue) - { - previous = (ulong)previousValue; - } - else - { - previous = ulong.MaxValue; - } - } - - ulong prime64 = (ulong)prime; - divisorData = MontgomeryDivisorData.FromModulus(prime64); - ulong order64 = CalculateGpu(gpu, prime64, previous, divisorData, config); - result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; - } - else - { - divisorData = MontgomeryDivisorData.Empty; - result = CalculateWideInternalGpu(gpu, prime, previousOrder, divisorData, config); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 CalculateWideInternalGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config) - { - if (prime <= UInt128.One) - { - return UInt128.One; - } - - if (prime == (UInt128)3UL) - { - return (UInt128)2UL; - } - - UInt128 phi = prime - UInt128.One; - PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); - UInt128 result; - if (phiFactors.Factors is null) - { - result = FinishStrictlyWideGpu(gpu, prime, divisorData); - } - else - { - result = RunHeuristicPipelineWideGpu(gpu, prime, previousOrder, divisorData, config, phi, phiFactors); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CollectFactorsGpu(PrimeOrderCalculatorAccelerator gpu, ref Span primeSlots, ref Span exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached) - { - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - - if (isPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); - continue; - } - - long currentTimestamp = Stopwatch.GetTimestamp(); - if (currentTimestamp > deadlineTimestamp) - { - pollardRhoDeadlineReached = true; - pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); - continue; - } - - // if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) - if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); - continue; - } - - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool EvaluateSpecialMaxCandidatesGpu( - PrimeOrderCalculatorAccelerator gpu, - ReadOnlySpan factors, - ulong phi, - ulong prime, - in MontgomeryDivisorData divisorData) - { - int factorCount = factors.Length; - gpu.EnsureUlongInputOutputCapacity(factorCount); - - int acceleratorIndex = gpu.AcceleratorIndex; - // GpuPrimeWorkLimiter.Acquire(); - ArrayView specialMaxFactorsView = gpu.InputView; - ArrayView specialMaxResultView = gpu.OutputUlongView2; - var kernelLauncher = gpu.SpecialMaxKernelLauncher; - - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - specialMaxFactorsView.SubView(0, factorCount).CopyFromCPU(stream, factors); - - kernelLauncher( - stream, - 1, - phi, - specialMaxFactorsView, - factorCount, - divisorData.Modulus, - specialMaxResultView); - - ulong result = 0UL; - specialMaxResultView.CopyToCPU(stream, ref result, 1); - stream.Synchronize(); - - AcceleratorStreamPool.Return(acceleratorIndex, stream); - // GpuPrimeWorkLimiter.Release(); - return result != 0; - } - - private static UInt128 ExponentLoweringWideGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors) - { - ArrayPool pool = ThreadStaticPools.FactorEntry128Pool; - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - // TODO: This should never trigger from production code - check - // if (length == 0) - // { - // return order; - // } - FactorEntry128[] tempArray = pool.Rent(length + 1); - - Span buffer = tempArray.AsSpan(0, length); - factorSpan.CopyTo(buffer); - - if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor)) - { - buffer[length] = new FactorEntry128(factors.Cofactor, 1); - length++; - } - - buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value)); - - // TODO: Investigating using stepper to replace modulo operation - for (int i = 0; i < length; i++) - { - UInt128 primeFactor = buffer[i].Value; - int exponent = buffer[i].Exponent; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) == UInt128.Zero) - { - UInt128 reduced = order / primeFactor; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - order = reduced; - continue; - } - } - - break; - } - } - - pool.Return(tempArray, clearArray: false); - return order; - } - - private static UInt128 FinishStrictlyWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData divisorData) - { - UInt128 phi = prime - UInt128.One; - Dictionary counts = new(capacity: 8); - FactorCompletelyWide(phi, counts); - if (counts.Count == 0) - { - return phi; - } - - List> entries = new(counts); - entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); - - UInt128 order = phi; - int entryCount = entries.Count; - for (int i = 0; i < entryCount; i++) - { - UInt128 primeFactor = entries[i].Key; - int exponent = entries[i].Value; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != UInt128.Zero) - { - UInt128 candidate = order / primeFactor; - if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) == UInt128.One) - { - order = candidate; - } - else - { - break; - } - } - } - } - - return order; - } - - private static UInt128 InitializeStartingOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData divisorData) - { - UInt128 order = phi; - UInt128 mod8 = prime & (UInt128)7UL; - if (mod8 == UInt128.One || mod8 == (UInt128)7UL) - { - UInt128 half = phi >> 1; - if (Pow2ModWideGpu(gpu, half, prime, divisorData) == UInt128.One) - { - order = half; - } - } - - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - // stackalloc is faster than pooling - var result = PartialFactorResult.Rent(divisorData); - ulong[] factors = result.Factors; - Span primeSlots = new (factors); - int[] exponents = result.Exponents; - Span exponentSlots = new (exponents); - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; - - bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out remaining); - - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); - - if (!gpuFactored) - { - PopulateSmallPrimeFactorsCpu( - value, - limit, - primeSlots, - exponentSlots, - out factorCount, - out remaining); - } - - bool limitReached; - if (remaining > 1UL) - { - if (config.PollardRhoMilliseconds > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - FixedCapacityStack compositeStack = ThreadStaticPools.RentUlongStack(4); - compositeStack.Push(remaining); - - CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); - - if (limitReached) - { - while (compositeStack.Count > 0) - { - pending.Add(new PartialFactorPendingEntry(compositeStack.Pop(), knownComposite: false)); - } - } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); - } - else - { - pending.Add(new PartialFactorPendingEntry(remaining, knownComposite: false)); - } - } - - ulong cofactor = 1UL; - bool cofactorContainsComposite = false; - int pendingCount = pending.Count; - int index = 0; - for (; index < pendingCount; index++) - { - PartialFactorPendingEntry entry = pending[index]; - ulong composite = entry.Value; - - if (entry.KnownComposite) - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - continue; - } - - if (!entry.HasKnownPrimality) - { - bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(composite); - - entry.WithPrimality(isPrime); - pending[index] = entry; - } - - if (entry.IsPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - } - } - - bool cofactorIsPrime; - if (cofactor <= 1UL) - { - cofactorIsPrime = false; - } - else if (cofactorContainsComposite) - { - cofactorIsPrime = false; - } - else - { - cofactorIsPrime = PrimeTesterByLastDigit.IsPrimeCpu(cofactor); - } - - // This will never happen in production code. We'll always get at least 1 factor - // if (factorCount == 0) - // { - // if (cofactor == value) - // { - // result = PartialFactorResult.Rent(null, value, false, 0); - // goto ReturnResult; - // } - - // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - // goto ReturnResult; - // } - - Array.Sort(factors, exponents, 0, factorCount); - - // Check if it was factored completely - result.Cofactor = cofactor; - result.FullyFactored = cofactor == 1UL; - result.Count = factorCount; - result.CofactorIsPrime = cofactorIsPrime; - - pending.Clear(); - ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 Pow2ModWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData) - { - if (modulus == UInt128.One) - { - return UInt128.Zero; - } - - if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue) - { - ulong prime64 = (ulong)modulus; - ulong exponent64 = (ulong)exponent; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent64, prime64, out ulong remainder, divisorData); - if (status == GpuPow2ModStatus.Success) - { - return remainder; - } - } - - GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, modulus, out UInt128 wideRemainder); - if (wideStatus == GpuPow2ModStatus.Success) - { - return wideRemainder; - } - - return exponent.Pow2MontgomeryModWindowed(modulus); - } - - private static UInt128 RunHeuristicPipelineWideGpu( - PrimeOrderCalculatorAccelerator gpu, - in UInt128 prime, - in UInt128? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - in UInt128 phi, - in PartialFactorResult128 phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxWideGpu(gpu, phi, prime, divisorData, phiFactors)) - { - return phi; - } - - UInt128 candidateOrder = InitializeStartingOrderWideGpu(gpu, prime, phi, divisorData); - candidateOrder = ExponentLoweringWideGpu(gpu, candidateOrder, prime, divisorData, phiFactors); - - if (TryConfirmOrderWideGpu(gpu, prime, candidateOrder, divisorData, config)) - { - return candidateOrder; - } - - if (config.StrictMode) - { - return FinishStrictlyWideGpu(gpu, prime, divisorData); - } - - if (TryHeuristicFinishWideGpu(gpu, prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) - { - return order; - } - - return candidateOrder; - } - - private static bool TryConfirmCandidateWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData) - { - PartialFactorResult128 factorization = PartialFactorWide(candidate, config); - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(factorization.Cofactor)) - { - return false; - } - - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - UInt128 primeFactor = span[i].Value; - UInt128 reduced = candidate; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if ((reduced % primeFactor) != UInt128.Zero) - { - break; - } - - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) - { - return false; - } - - powUsed++; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - } - - return true; - } - - private static bool TryConfirmOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config) - { - if (order == UInt128.Zero) - { - return false; - } - - if (Pow2ModWideGpu(gpu, order, prime, divisorData) != UInt128.One) - { - return false; - } - - PartialFactorResult128 factorization = PartialFactorWide(order, config); - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(factorization.Cofactor)) - { - return false; - } - - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - UInt128 primeFactor = span[i].Value; - UInt128 reduced = order; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) - { - break; - } - - reduced /= primeFactor; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - } - - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryHeuristicFinishGpu( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - PartialFactorResult phiFactors, - PartialFactorResult? cachedOrderFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - cachedOrderFactors?.Dispose(); - return false; - } - - // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, divisorData, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime) - { - return false; - } - - orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); - candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - ThreadStaticPools.ReturnUlongList(candidates); - return false; - } - - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - bool allowGpuBatch = true; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool powStepperInitialized = false; - - int index = 0; - Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; - ArrayPool pool = ThreadStaticPools.UlongPool; - while (index < candidateCount && powUsed < powBudget) - { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize)); - if (batchSize <= 0) - { - break; - } - - ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); - ulong[]? gpuPool = null; - Span pooledGpuRemainders = default; - bool gpuSuccess = false; - bool gpuStackRemainders = false; - GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; - - if (allowGpuBatch) - { - if (batchSize <= PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount) - { - Span localRemainders = stackGpuRemainders[..batchSize]; - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, localRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - gpuSuccess = true; - gpuStackRemainders = true; - } - } - else - { - gpuPool = pool.Rent(batchSize); - Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, pooledRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - pooledGpuRemainders = pooledRemainders; - gpuSuccess = true; - } - else - { - pool.Return(gpuPool, clearArray: false); - gpuPool = null; - } - } - - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) - { - allowGpuBatch = false; - } - } - - int startIndex = 0; - - if (!gpuSuccess && !powStepperInitialized) - { - ulong candidate = batch[0]; - powUsed++; - - bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); - powStepperInitialized = true; - if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - startIndex = 1; - } - - for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) - { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne; - if (gpuSuccess) - { - ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; - equalsOne = remainderValue == 1UL; - } - else - { - equalsOne = powStepper.ComputeNextIsUnity(candidate); - } - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - index += batchSize; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - // DebugLog("No candidate confirmed"); - return false; - } - finally - { - orderFactors.Dispose(); - } - } - - private static bool TryHeuristicFinishWideGpu( - PrimeOrderCalculatorAccelerator gpu, - in UInt128 prime, - in UInt128 order, - in UInt128? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - out UInt128 result) - { - result = UInt128.Zero; - if (order <= UInt128.One) - { - return false; - } - - PartialFactorResult128 orderFactors = PartialFactorWide(order, config); - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(orderFactors.Cofactor)) - { - return false; - } - - orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = new(capacity); - FactorEntry128[] factorArray = orderFactors.Factors!; - BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); - if (candidates.Count == 0) - { - return false; - } - - SortCandidatesWide(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; - int powUsed = 0; - int candidateCount = candidates.Count; - - for (int i = 0; i < candidateCount; i++) - { - if (powUsed >= powBudget && powBudget > 0) - { - break; - } - - UInt128 candidate = candidates[i]; - powUsed++; - - if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) != UInt128.One) - { - continue; - } - - if (!TryConfirmCandidateWideGpu(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData)) - { - continue; - } - - result = candidate; - return true; - } - - return false; - } - - private static bool TryPollardRhoGpu(PrimeOrderCalculatorAccelerator gpu, ulong n, out ulong factor) - { - Span randomStateSpan = stackalloc ulong[1]; - randomStateSpan[0] = ThreadStaticDeterministicRandomGpu.Exclusive.State; - - Span factoredSpan = stackalloc byte[1]; - Span factorSpan = stackalloc ulong[1]; - - int acceleratorIndex = gpu.AcceleratorIndex; - var kernel = gpu.PollardRhoKernel; - - AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); - gpu.InputView.CopyFromCPU(stream, randomStateSpan); - - kernel.Launch(stream, 1, n, 1, gpu.InputView, gpu.OutputByteView, gpu.OutputUlongView); - - gpu.OutputByteView.CopyToCPU(stream, factoredSpan); - gpu.OutputUlongView.CopyToCPU(stream, factorSpan); - gpu.InputView.CopyToCPU(stream, randomStateSpan); - stream.Synchronize(); - - AcceleratorStreamPool.Return(acceleratorIndex, stream); - ThreadStaticDeterministicRandomGpu.Exclusive.SetState(randomStateSpan[0]); - - bool factored = factoredSpan[0] != 0; - factor = factored ? factorSpan[0] : 0UL; - return factored; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryPopulateSmallPrimeFactorsGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, uint limit, Dictionary counts, out int factorCount, out ulong remaining) - { - var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); - var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); - Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - remaining = value; - - gpu.EnsureSmallPrimeFactorSlotsCapacity(GpuSmallPrimeFactorSlots); - int acceleratorIndex = gpu.AcceleratorIndex; - var accelerator = gpu.Accelerator; - - var kernelLauncher = gpu.SmallPrimeFactorKernelLauncher; - ArrayView1D smallPrimeFactorCountSlotView = gpu.OutputIntView2; - ArrayView1D smallPrimeFactorRemainingSlotView = gpu.InputView; - ArrayView1D smallPrimeFactorPrimeSlotsView = gpu.OutputUlongView; - ArrayView1D smallPrimeFactorExponentSlotsView = gpu.OutputIntView; - ArrayView1D smallPrimeFactorsSquaresView = gpu.SmallPrimeFactorSquares; - ArrayView1D smallPrimeFactorsPrimesView = gpu.SmallPrimeFactorPrimes; - - // GpuPrimeWorkLimiter.Acquire(); - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - kernelLauncher( - stream, - 1, - value, - limit, - smallPrimeFactorsPrimesView, - smallPrimeFactorsSquaresView, - (int)smallPrimeFactorsPrimesView.Length, - smallPrimeFactorPrimeSlotsView, - smallPrimeFactorExponentSlotsView, - smallPrimeFactorCountSlotView, - smallPrimeFactorRemainingSlotView); - - factorCount = 0; - smallPrimeFactorCountSlotView.CopyToCPU(stream, ref factorCount, 1); - stream.Synchronize(); - - factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); - - if (factorCount > 0) - { - primeBuffer = primeBuffer[..factorCount]; - smallPrimeFactorPrimeSlotsView.CopyToCPU(stream, primeBuffer); - - exponentBuffer = exponentBuffer[..factorCount]; - smallPrimeFactorExponentSlotsView.CopyToCPU(stream, exponentBuffer); - } - - smallPrimeFactorRemainingSlotView.CopyToCPU(stream, ref remaining, 1); - stream.Synchronize(); - - AcceleratorStreamPool.Return(acceleratorIndex, stream); - // GpuPrimeWorkLimiter.Release(); - - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - counts.Add(primeValue, exponent); - } - - return true; - } - - private static bool TrySpecialMaxWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors) - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - for (int i = 0; i < length; i++) - { - UInt128 factor = factorSpan[i].Value; - UInt128 reduced = phi / factor; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - - return true; - } -} - diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs deleted file mode 100644 index e0b936e..0000000 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs +++ /dev/null @@ -1,651 +0,0 @@ -using System.Buffers; -using System.Diagnostics; -using System.Runtime.CompilerServices; -using System.Runtime.InteropServices; -using ILGPU; -using PerfectNumbers.Core.Gpu; -using PerfectNumbers.Core.Gpu.Accelerators; - -namespace PerfectNumbers.Core; - -internal static partial class PrimeOrderCalculator -{ - public static ulong CalculateHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config) - { - // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. - if (prime <= 3UL) - { - return prime == 3UL ? 2UL : 1UL; - } - - ulong phi = prime - 1UL; - - PartialFactorResult phiFactors = PartialFactorHybrid(gpu, phi, divisorData, config); - - ulong result; - if (phiFactors.Factors is null) - { - result = CalculateByFactorizationCpu(prime, divisorData); - - phiFactors.Dispose(); - return result; - } - - result = RunHeuristicPipelineHybrid(gpu, prime, previousOrder, config, divisorData, phi, phiFactors); - - phiFactors.Dispose(); - return result; - } - - public static UInt128 CalculateHybrid( - PrimeOrderCalculatorAccelerator gpu, - in UInt128 prime, - in UInt128? previousOrder, - in PrimeOrderCalculatorConfig config) - { - MontgomeryDivisorData divisorData; - UInt128 result; - if (prime <= ulong.MaxValue) - { - ulong? previous = null; - if (previousOrder.HasValue) - { - UInt128 previousValue = previousOrder.Value; - if (previousValue <= ulong.MaxValue) - { - previous = (ulong)previousValue; - } - else - { - previous = ulong.MaxValue; - } - } - - ulong prime64 = (ulong)prime; - divisorData = MontgomeryDivisorData.FromModulus(prime64); - ulong order64 = CalculateHybrid(gpu, prime64, previous, divisorData, config); - result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; - } - else - { - divisorData = MontgomeryDivisorData.Empty; - result = CalculateWideInternalCpu(prime, previousOrder, divisorData, config); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span primeSlots, Span exponentSlots, ref int factorCount, List pending, FixedCapacityStack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) - { - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - - if (isPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); - continue; - } - - pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; - if (pollardRhoDeadlineReached) - { - pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); - continue; - } - - if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); - continue; - } - - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - - pollardRhoDeadlineReached = false; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong InitializeStartingOrderHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong phi, in MontgomeryDivisorData divisorData) - { - ulong order = phi; - if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) - { - ulong half = phi >> 1; - if (Pow2EqualsOneHybrid(gpu, half, prime, divisorData)) - { - order = half; - } - } - - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - // stackalloc is faster than pooling - - var result = PartialFactorResult.Rent(divisorData); - ulong[] factors = result.Factors; - Span primeSlots = new(factors); - int[] exponents = result.Exponents; - Span exponentSlots = new(exponents); - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); - - PopulateSmallPrimeFactorsCpu( - value, - limit, - primeSlots, - exponentSlots, - out int factorCount, - out ulong remaining); - - bool limitReached; - if (remaining > 1UL) - { - if (config.PollardRhoMilliseconds > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - FixedCapacityStack compositeStack = ThreadStaticPools.RentUlongStack(4); - compositeStack.Push(remaining); - - CollectFactorsHybrid(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); - - if (limitReached) - { - while (compositeStack.Count > 0) - { - pending.Add(new PartialFactorPendingEntry(compositeStack.Pop(), knownComposite: false)); - } - } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); - } - else - { - pending.Add(new PartialFactorPendingEntry(remaining, knownComposite: false)); - } - } - - ulong cofactor = 1UL; - bool cofactorContainsComposite = false; - int pendingCount = pending.Count; - int index = 0; - for (; index < pendingCount; index++) - { - PartialFactorPendingEntry entry = pending[index]; - ulong composite = entry.Value; - - if (entry.KnownComposite) - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - continue; - } - - if (!entry.HasKnownPrimality) - { - bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - - entry.WithPrimality(isPrime); - pending[index] = entry; - } - - if (entry.IsPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - } - } - - bool cofactorIsPrime; - if (cofactor <= 1UL) - { - cofactorIsPrime = false; - } - else if (cofactorContainsComposite) - { - cofactorIsPrime = false; - } - else - { - cofactorIsPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(cofactor) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor); - } - - // This will never happen in production code. We'll always get at least 1 factor - // if (factorCount == 0) - // { - // if (cofactor == value) - // { - // result = PartialFactorResult.Rent(null, value, false, 0); - // goto ReturnResult; - // } - - // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - // goto ReturnResult; - // } - - Array.Sort(factors, exponents, 0, factorCount); - - result.Cofactor = cofactor; - result.FullyFactored = cofactor == 1UL; - result.Count = factorCount; - result.CofactorIsPrime = cofactorIsPrime; - - pending.Clear(); - ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static bool Pow2EqualsOneHybrid(PrimeOrderCalculatorAccelerator gpu, ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) - { - // if (exponent >= PerfectNumberConstants.MaxQForDivisorCycles || prime >= PerfectNumberConstants.MaxQForDivisorCycles) - // { - // // Atomic.Add(ref _pow2GpuHits, 1); - // // Console.WriteLine($"pow2 GPU Hits {_pow2GpuHits}"); - // ulong remainder = exponent.Pow2MontgomeryModWindowedConvertGpu(gpu, divisorData); - // // GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); - // // if (status == GpuPow2ModStatus.Success) - // // { - // return remainder == 1UL; - // // } - // } - - // Atomic.Add(ref _pow2CpuHits, 1); - // Console.WriteLine($"pow2 CPU Hits {_pow2CpuHits}"); - return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong RunHeuristicPipelineHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong? previousOrder, - in PrimeOrderCalculatorConfig config, - in MontgomeryDivisorData divisorData, - ulong phi, - PartialFactorResult phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, prime, phiFactors, divisorData)) - { - return phi; - } - - ulong candidateOrder = InitializeStartingOrderHybrid(gpu, prime, phi, divisorData); - candidateOrder = ExponentLoweringCpu(candidateOrder, phiFactors, divisorData); - - if (TryConfirmOrderHybrid(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) - { - return candidateOrder; - } - - if (config.StrictMode) - { - orderFactors?.Dispose(); - return CalculateByFactorizationCpu(prime, divisorData); - } - - if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) - { - return order; - } - - return candidateOrder; - } - - private static int _cpuCount; - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool RunOnCpu() - { - int gpuRatio = EnvironmentConfiguration.GpuRatio; - if (gpuRatio == 0) - { - return false; - } - - if (gpuRatio == 1) - { - while (true) - { - if (Interlocked.CompareExchange(ref _cpuCount, 0, 1) == 0) - { - return true; - } - - if (Interlocked.CompareExchange(ref _cpuCount, 1, 0) == 1) - { - return false; - } - - Thread.Yield(); - } - } - - int cpuCount = Interlocked.Increment(ref _cpuCount); - if (cpuCount == gpuRatio) - { - Interlocked.Add(ref _cpuCount, -gpuRatio); - } - else if (cpuCount > gpuRatio) - { - cpuCount -= gpuRatio; - } - - return cpuCount != gpuRatio; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) - { - PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, divisorData, config); - - if (factorization.Factors is null) - { - factorization.Dispose(); - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) - { - factorization.Dispose(); - return false; - } - - factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan factorsSpan = factorization.Factors; - ReadOnlySpan exponentsSpan = factorization.Exponents; - int length = factorization.Count; - - const int StackExponentCapacity = 32; - const int ExponentHardLimit = 256; - - Span stackBuffer = stackalloc ulong[StackExponentCapacity]; - ulong[]? heapCandidateArray = null; - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool violates = false; - - for (int i = 0; i < length; i++) - { - ulong primeFactor = factorsSpan[i]; - int exponent = exponentsSpan[i]; - - // Factorization entries for candidate confirmation are always positive here. - // if (exponent <= 0) - // { - // continue; - // } - - if (exponent > ExponentHardLimit) - { - throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); - } - - if (exponent <= StackExponentCapacity) - { - if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) - { - violates = true; - break; - } - - continue; - } - - heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit); - if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) - { - violates = true; - break; - } - } - - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - if (heapCandidateArray is not null) - { - ThreadStaticPools.UlongPool.Return(heapCandidateArray); - } - factorization.Dispose(); - - return !violates; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmOrderHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong order, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - out PartialFactorResult? reusableFactorization) - { - reusableFactorization = null; - - if (order == 0UL) - { - return false; - } - - if (!Pow2EqualsOneHybrid(gpu, order, prime, divisorData)) - { - return false; - } - - PartialFactorResult factorization = PartialFactorHybrid(gpu, order, divisorData, config); - if (factorization.Factors is null) - { - factorization.Dispose(); - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) - { - factorization.Dispose(); - return false; - } - - factorization.WithAdditionalPrime(factorization.Cofactor); - } - - if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) - { - reusableFactorization = factorization; - return false; - } - - factorization.Dispose(); - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryHeuristicFinishHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - PartialFactorResult phiFactors, - PartialFactorResult? cachedOrderFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - cachedOrderFactors?.Dispose(); - return false; - } - - // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, divisorData, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime) - { - return false; - } - - orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); - candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - ThreadStaticPools.ReturnUlongList(candidates); - return false; - } - - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool powStepperInitialized = false; - - int index = 0; - Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; - ArrayPool pool = ThreadStaticPools.UlongPool; - while (index < candidateCount && powUsed < powBudget) - { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize)); - if (batchSize <= 0) - { - break; - } - - ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); - int startIndex = 0; - - if (!powStepperInitialized) - { - ulong candidate = batch[0]; - powUsed++; - - bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); - powStepperInitialized = true; - if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - startIndex = 1; - } - - for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) - { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne = powStepper.ComputeNextIsUnity(candidate); - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - index += batchSize; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - - return false; - } - finally - { - orderFactors.Dispose(); - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - int length = factors.Count; - // The phi factorization on the scanner path always yields at least one entry for phi >= 2, - // preventing the zero-length case from occurring outside synthetic tests. - // if (length == 0) - // { - // return true; - // } - - ReadOnlySpan factorSpan = new(factors.Factors, 0, length); - Span stackBuffer = stackalloc ulong[length]; - - if (length <= 8) - { - return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); - } - - return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData); - } - -} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs deleted file mode 100644 index 3681aec..0000000 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ /dev/null @@ -1,631 +0,0 @@ -using System.Buffers; -using System.Diagnostics; -using System.Numerics; - -namespace PerfectNumbers.Core; - -internal static partial class PrimeOrderCalculator -{ - private static void SortCandidatesWide(UInt128 prime, in UInt128? previousOrder, List candidates) - { - UInt128 previous = previousOrder ?? UInt128.Zero; - int previousGroup = previousOrder.HasValue ? GetGroupWide(previousOrder.Value, prime) : 1; - bool hasPrevious = previousOrder.HasValue; - - candidates.Sort((x, y) => - { - CandidateKey128 keyX = BuildKeyWide(x, prime, previous, previousGroup, hasPrevious); - CandidateKey128 keyY = BuildKeyWide(y, prime, previous, previousGroup, hasPrevious); - - int primary = keyX.Primary.CompareTo(keyY.Primary); - if (primary != 0) - { - return primary; - } - - int secondary = CompareComponents(keyX.SecondaryDescending, keyX.Secondary, keyY.SecondaryDescending, keyY.Secondary); - if (secondary != 0) - { - return secondary; - } - - return CompareComponents(keyX.TertiaryDescending, keyX.Tertiary, keyY.TertiaryDescending, keyY.Tertiary); - }); - } - - private static CandidateKey128 BuildKeyWide( - in UInt128 value, - in UInt128 prime, - in UInt128 previous, - int previousGroup, - bool hasPrevious) - { - int group = GetGroupWide(value, prime); - if (group == 0) - { - return new CandidateKey128(int.MaxValue, false, UInt128.Zero, false, UInt128.Zero); - } - - bool isGe = !hasPrevious || value >= previous; - int primary = ComputePrimary(group, isGe, previousGroup); - - if (group == 3) - { - return new CandidateKey128(primary, true, value, true, value); - } - - UInt128 distance = hasPrevious ? (value > previous ? value - previous : previous - value) : value; - return new CandidateKey128(primary, false, distance, false, value); - } - - private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY) - { - if (descendingX == descendingY) - { - if (valueX == valueY) - { - return 0; - } - - if (descendingX) - { - return valueX > valueY ? -1 : 1; - } - - return valueX > valueY ? 1 : -1; - } - - if (valueX == valueY) - { - return 0; - } - - return valueX > valueY ? 1 : -1; - } - - private static int GetGroupWide(in UInt128 value, in UInt128 prime) - { - UInt128 threshold1 = prime >> 3; - if (value <= threshold1) - { - return 1; - } - - UInt128 threshold2 = prime >> 2; - if (value <= threshold2) - { - return 2; - } - - UInt128 threshold3 = (UInt128)(((BigInteger)prime * 3) >> 3); - if (value <= threshold3) - { - return 3; - } - - return 0; - } - - private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List candidates, int limit) - { - if (count == 0) - { - return; - } - - Span buffer = factors.AsSpan(0, count); - buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - BuildCandidatesRecursiveWide(order, buffer, 0, UInt128.One, candidates, limit); - } - - private static void BuildCandidatesRecursiveWide( - in UInt128 order, - in ReadOnlySpan factors, - int index, - in UInt128 divisorProduct, - List candidates, - int limit) - { - if (candidates.Count >= limit) - { - return; - } - - if (index >= factors.Length) - { - if (divisorProduct == UInt128.One || divisorProduct == order) - { - return; - } - - UInt128 candidate = order / divisorProduct; - if (candidate > UInt128.One && candidate < order) - { - candidates.Add(candidate); - } - - return; - } - - FactorEntry128 factor = factors[index]; - UInt128 primeFactor = factor.Value; - UInt128 contribution = UInt128.One; - for (int exponent = 0; exponent <= factor.Exponent; exponent++) - { - UInt128 nextDivisor = divisorProduct * contribution; - if (nextDivisor > order) - { - break; - } - - BuildCandidatesRecursiveWide(order, factors, index + 1, nextDivisor, candidates, limit); - if (candidates.Count >= limit) - { - return; - } - - if (exponent == factor.Exponent) - { - break; - } - - if (contribution > order / primeFactor) - { - break; - } - - contribution *= primeFactor; - } - } - - private static PartialFactorResult128 PartialFactorWide(in UInt128 value, in PrimeOrderCalculatorConfig config) - { - if (value <= UInt128.One) - { - return PartialFactorResult128.Empty; - } - - Dictionary counts = new(capacity: 8); - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - UInt128 remaining; - - // if (IsGpuPow2Allowed && value <= ulong.MaxValue) - // { - // Dictionary narrowCounts = new(capacity: 8); - // if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining)) - // { - // foreach (KeyValuePair entry in narrowCounts) - // { - // counts[(UInt128)entry.Key] = entry.Value; - // } - - // remaining = (UInt128)narrowRemaining; - // } - // else - // { - // remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - // } - // } - // else - // { - remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - // } - - List pending = new(); - if (remaining > UInt128.One) - { - pending.Add(remaining); - } - - if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - Stack stack = new(); - stack.Push(remaining); - pending.Clear(); - - long timestamp = 0L; // reused for deadline checks. - while (stack.Count > 0) - { - UInt128 composite = stack.Pop(); - if (composite <= UInt128.One) - { - continue; - } - - if (IsPrimeWide(composite)) - { - AddFactor(counts, composite, 1); - continue; - } - - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - pending.Add(composite); - continue; - } - - if (!TryPollardRhoWide(composite, deadlineTimestamp, out UInt128 factor)) - { - pending.Add(composite); - continue; - } - - UInt128 quotient = composite / factor; - stack.Push(factor); - stack.Push(quotient); - } - } - - UInt128 cofactor = UInt128.One; - int pendingCount = pending.Count; - for (int i = 0; i < pendingCount; i++) - { - UInt128 composite = pending[i]; - if (IsPrimeWide(composite)) - { - AddFactor(counts, composite, 1); - } - else - { - BigInteger product = (BigInteger)cofactor * (BigInteger)composite; - cofactor = (UInt128)product; - } - } - - if (counts.Count == 0 && cofactor == value) - { - return new PartialFactorResult128(null, value, false, 0); - } - - ArrayPool? pool = null; - FactorEntry128[]? rented = counts.Count > 0 ? (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null; - int index = 0; - if (rented is not null) - { - foreach (KeyValuePair entry in counts) - { - rented[index] = new FactorEntry128(entry.Key, entry.Value); - index++; - } - } - - FactorEntry128[]? resultArray = null; - if (pool is not null) - { - Span span = rented.AsSpan(0, index); - span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - resultArray = new FactorEntry128[index]; - span.CopyTo(resultArray); - pool.Return(rented!, clearArray: false); - } - - bool fullyFactored = cofactor == UInt128.One; - return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index); - } - - private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary 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(in UInt128 n, long deadlineTimestamp, out UInt128 factor) - { - factor = UInt128.Zero; - if ((n & UInt128.One) == UInt128.Zero) - { - factor = (UInt128)2UL; - return true; - } - - UInt128 c = UInt128.One; - UInt128 x = (UInt128)2UL; - UInt128 y = x; - long timestamp = 0L; // reused for deadline checks. - - while (true) - { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } - - x = AdvancePolynomialWide(x, c, n); - y = AdvancePolynomialWide(y, c, n); - y = AdvancePolynomialWide(y, c, n); - - UInt128 diff = x > y ? x - y : y - x; - UInt128 d = BinaryGcdWide(diff, n); - if (d == n) - { - c += UInt128.One; - x = (UInt128)2UL; - y = x; - continue; - } - - if (d > UInt128.One) - { - factor = d; - return true; - } - } - } - - private static UInt128 AdvancePolynomialWide(in UInt128 x, in UInt128 c, in UInt128 modulus) - { - BigInteger value = (BigInteger)x; - value = (value * value + (BigInteger)c) % (BigInteger)modulus; - return (UInt128)value; - } - - private static UInt128 BinaryGcdWide(UInt128 a, UInt128 b) - { - if (a == UInt128.Zero) - { - return b; - } - - if (b == UInt128.Zero) - { - return a; - } - - int shift = TrailingZeroCountWide(a | b); - a >>= TrailingZeroCountWide(a); - - while (true) - { - b >>= TrailingZeroCountWide(b); - if (a > b) - { - (a, b) = (b, a); - } - - b -= a; - if (b == UInt128.Zero) - { - return a << shift; - } - } - } - - private static int TrailingZeroCountWide(in UInt128 value) - { - if (value == UInt128.Zero) - { - return 128; - } - - ulong low = (ulong)value; - if (low != 0UL) - { - return BitOperations.TrailingZeroCount(low); - } - - ulong high = (ulong)(value >> 64); - return 64 + BitOperations.TrailingZeroCount(high); - } - - private static void FactorCompletelyWide(in UInt128 value, Dictionary counts) - { - if (value <= UInt128.One) - { - return; - } - - if (IsPrimeWide(value)) - { - AddFactor(counts, value, 1); - return; - } - - UInt128 factor = PollardRhoWide(value); - FactorCompletelyWide(factor, counts); - FactorCompletelyWide(value / factor, counts); - } - - private static UInt128 PollardRhoWide(in UInt128 n) - { - if ((n & UInt128.One) == UInt128.Zero) - { - return (UInt128)2UL; - } - - UInt128 c = UInt128.One; - while (true) - { - UInt128 x = (UInt128)2UL; - UInt128 y = x; - UInt128 d = UInt128.One; - - while (d == UInt128.One) - { - x = AdvancePolynomialWide(x, c, n); - y = AdvancePolynomialWide(y, c, n); - y = AdvancePolynomialWide(y, c, n); - - UInt128 diff = x > y ? x - y : y - x; - d = BinaryGcdWide(diff, n); - } - - if (d != n) - { - return d; - } - - c += UInt128.One; - } - } - - private static bool IsPrimeWide(in UInt128 value) - { - if (value <= ulong.MaxValue) - { - // return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None); - return PrimeTester.IsPrimeCpu((ulong)value); - } - - if ((value & UInt128.One) == UInt128.Zero) - { - return false; - } - - uint[] smallPrimes = PrimesGenerator.SmallPrimes; - int length = smallPrimes.Length; - for (int i = 0; i < length; i++) - { - uint prime = smallPrimes[i]; - UInt128 primeValue = prime; - UInt128 primeSquare = primeValue * primeValue; - if (primeSquare > value) - { - break; - } - - if ((value % primeValue) == UInt128.Zero) - { - return false; - } - } - - BigInteger n = (BigInteger)value; - BigInteger d = n - BigInteger.One; - int s = 0; - while ((d & 1) == 0) - { - d >>= 1; - s++; - } - - int[] bases = new[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 }; - for (int i = 0; i < bases.Length; i++) - { - int baseValue = bases[i]; - if (baseValue <= 1) - { - continue; - } - - BigInteger a = baseValue; - if (a >= n) - { - continue; - } - - BigInteger x = BigInteger.ModPow(a, d, n); - if (x == BigInteger.One || x == n - BigInteger.One) - { - continue; - } - - bool witnessFound = true; - for (int r = 1; r < s; r++) - { - x = BigInteger.ModPow(x, 2, n); - if (x == n - BigInteger.One) - { - witnessFound = false; - break; - } - } - - if (witnessFound) - { - return false; - } - } - - return true; - } - - private static void AddFactor(Dictionary counts, in UInt128 prime, int exponent) - { - if (counts.TryGetValue(prime, out int existing)) - { - counts[prime] = existing + exponent; - } - else - { - counts[prime] = exponent; - } - } - - private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary) - { - public int Primary { get; } = primary; - - public bool SecondaryDescending { get; } = secondaryDescending; - - public UInt128 Secondary { get; } = secondary; - - public bool TertiaryDescending { get; } = tertiaryDescending; - - public UInt128 Tertiary { get; } = tertiary; - } - - private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count) - { - public FactorEntry128[]? Factors { get; } = factors; - - public UInt128 Cofactor { get; } = cofactor; - - public bool FullyFactored { get; } = fullyFactored; - - public int Count { get; } = count; - - public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0); - - public PartialFactorResult128 WithAdditionalPrime(UInt128 prime) - { - if (Factors is null) - { - FactorEntry128[] local = new FactorEntry128[1]; - local[0] = new FactorEntry128(prime, 1); - return new PartialFactorResult128(local, UInt128.One, true, 1); - } - - FactorEntry128[] extended = new FactorEntry128[Count + 1]; - Array.Copy(Factors, extended, Count); - extended[Count] = new FactorEntry128(prime, 1); - return new PartialFactorResult128(extended, UInt128.One, true, Count + 1); - } - } -} diff --git a/PerfectNumbers.Core/PrimeOrderCalculatorConfig.cs b/PerfectNumbers.Core/PrimeOrderCalculatorConfig.cs new file mode 100644 index 0000000..c363f88 --- /dev/null +++ b/PerfectNumbers.Core/PrimeOrderCalculatorConfig.cs @@ -0,0 +1,34 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +public readonly struct PrimeOrderCalculatorConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, bool strictMode, int maxSeconds = 0) +{ + private static PrimeOrderCalculatorConfig _heuristicDefault = CreateHeuristicDefault(UnboundedTaskScheduler.ConfiguredThreadCount); + + public static PrimeOrderCalculatorConfig HeuristicDefault + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get => _heuristicDefault; + } + + public static readonly PrimeOrderCalculatorConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, true); + + public static void ConfigureHeuristicDefault(int threadCount) + { + _heuristicDefault = CreateHeuristicDefault(threadCount); + } + + public static PrimeOrderCalculatorConfig CreateHeuristicDefault(int threadCount) + { + int normalizedThreadCount = threadCount < 1 ? 1 : threadCount; + return new(smallFactorLimit: 100_000, pollardRhoMilliseconds: normalizedThreadCount * 24, maxPowChecks: 24, false); + } + + public readonly uint SmallFactorLimit = smallFactorLimit; + public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; + public readonly int MaxPowChecks = maxPowChecks; + public readonly int MaxPowChecksCapacity = maxPowChecks << 2; + public readonly bool StrictMode = strictMode; + public readonly int MaxSeconds = maxSeconds; +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculatorTemplate.cs b/PerfectNumbers.Core/PrimeOrderCalculatorTemplate.cs new file mode 100644 index 0000000..f8bc931 --- /dev/null +++ b/PerfectNumbers.Core/PrimeOrderCalculatorTemplate.cs @@ -0,0 +1,3273 @@ +using System.Buffers; +using System.Collections.Generic; +using System.Collections.Concurrent; +using System.Diagnostics; +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; +using System.Threading; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core.Cpu; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core; + +[DeviceDependentTemplate(typeof(ComputationDevice))] +public static partial class PrimeOrderCalculatorTemplate +{ + // TODO: Remove branching on the CPU / GPU paths + // TODO: Split big non-static functions into smaller / extract static code to limit JIT time + + private static readonly ConcurrentDictionary _dictionaryStats = new(20_480, 64); + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static ulong BinaryGcd(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + a >>= BitOperations.TrailingZeroCount(a); + + while (true) + { + b >>= BitOperations.TrailingZeroCount(b); + if (a > b) + { + (a, b) = (b, a); + } + + b -= a; + if (b == 0UL) + { + return a << shift; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static ulong PollardRhoStrict(ulong n) + { + if ((n & 1UL) == 0UL) + { + return 2UL; + } + + while (true) + { + ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomialTwice(y, c, n); + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, n); + } + + if (d != n) + { + return d; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static int ComputePrimary(int group, bool isGe, int previousGroup) + { + int groupOffset = group switch + { + 1 => 0, + 2 => 2, + 3 => 4, + _ => 6, + }; + + if (group == previousGroup) + { + if (group == 3) + { + return groupOffset + (isGe ? 0 : 3); + } + + return groupOffset + (isGe ? 0 : 1); + } + + return groupOffset + (isGe ? 0 : 1); + } + + private static void PrintStats(int count) + { + _dictionaryStats.AddOrUpdate(count, 1UL, static (_, existing) => existing + 1UL); + Console.WriteLine("Current stats:"); + foreach (KeyValuePair entry in _dictionaryStats) + { + Console.WriteLine($"Count {entry.Key} = {entry.Value}"); + } + } + + private static CandidateKey128 BuildKeyWide( + in UInt128 value, + in UInt128 prime, + in UInt128 previous, + int previousGroup, + bool hasPrevious) + { + int group = GetGroupWide(value, prime); + if (group == 0) + { + return new CandidateKey128(int.MaxValue, false, UInt128.Zero, false, UInt128.Zero); + } + + bool isGe = !hasPrevious || value >= previous; + int primary = ComputePrimary(group, isGe, previousGroup); + + if (group == 3) + { + return new CandidateKey128(primary, true, value, true, value); + } + + UInt128 distance = hasPrevious ? (value > previous ? value - previous : previous - value) : value; + return new CandidateKey128(primary, false, distance, false, value); + } + + private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY) + { + if (descendingX == descendingY) + { + if (valueX == valueY) + { + return 0; + } + + if (descendingX) + { + return valueX > valueY ? -1 : 1; + } + + return valueX > valueY ? 1 : -1; + } + + if (valueX == valueY) + { + return 0; + } + + return valueX > valueY ? 1 : -1; + } + + private static int GetGroupWide(in UInt128 value, in UInt128 prime) + { + UInt128 threshold1 = prime >> 3; + if (value <= threshold1) + { + return 1; + } + + UInt128 threshold2 = prime >> 2; + if (value <= threshold2) + { + return 2; + } + + UInt128 threshold3 = (UInt128)(((BigInteger)prime * 3) >> 3); + if (value <= threshold3) + { + return 3; + } + + return 0; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static void BuildCandidates(ulong order, ulong[] factors, int[] exponents, int count, List candidates, int limit) + { + if (count == 0) + { + return; + } + + Array.Sort(factors, exponents, 0, count); + BuildCandidatesRecursive(order, factors, exponents, 0, 1UL, candidates, limit); + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void BuildCandidatesRecursive(ulong order, in ulong[] factors, in int[] exponents, int index, ulong divisorProduct, List candidates, int limit) + { + if (candidates.Count >= limit) + { + return; + } + + if (index >= factors.Length) + { + if (divisorProduct == 1UL || divisorProduct == order) + { + return; + } + + ulong candidate = order / divisorProduct; + if (candidate > 1UL && candidate < order) + { + candidates.Add(candidate); + } + + return; + } + + // FactorEntry factor = factors[index]; + // int factorExponent = factor.Exponent; + // ulong primeFactor = factor.Value; + + ulong factor = factors[index]; + int factorExponent = exponents[index]; + + ulong contribution = 1UL; + ulong contributionLimit = factorExponent == 0 ? 0UL : order / factor; + for (int exponent = 0; exponent <= factorExponent; exponent++) + { + ulong nextDivisor = divisorProduct * contribution; + if (nextDivisor > order) + { + break; + } + + BuildCandidatesRecursive(order, factors, exponents, index + 1, nextDivisor, candidates, limit); + if (candidates.Count >= limit) + { + return; + } + + if (exponent == factorExponent) + { + break; + } + + if (contribution > contributionLimit) + { + break; + } + + contribution *= factor; + } + } + + private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List candidates, int limit) + { + if (count == 0) + { + return; + } + + Span buffer = factors.AsSpan(0, count); + buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + BuildCandidatesRecursiveWide(order, buffer, 0, UInt128.One, candidates, limit); + } + + private static void BuildCandidatesRecursiveWide( + in UInt128 order, + in ReadOnlySpan factors, + int index, + in UInt128 divisorProduct, + List candidates, + int limit) + { + if (candidates.Count >= limit) + { + return; + } + + if (index >= factors.Length) + { + if (divisorProduct == UInt128.One || divisorProduct == order) + { + return; + } + + UInt128 candidate = order / divisorProduct; + if (candidate > UInt128.One && candidate < order) + { + candidates.Add(candidate); + } + + return; + } + + FactorEntry128 factor = factors[index]; + UInt128 primeFactor = factor.Value; + UInt128 contribution = UInt128.One; + for (int exponent = 0; exponent <= factor.Exponent; exponent++) + { + UInt128 nextDivisor = divisorProduct * contribution; + if (nextDivisor > order) + { + break; + } + + BuildCandidatesRecursiveWide(order, factors, index + 1, nextDivisor, candidates, limit); + if (candidates.Count >= limit) + { + return; + } + + if (exponent == factor.Exponent) + { + break; + } + + if (contribution > order / primeFactor) + { + break; + } + + contribution *= primeFactor; + } + } + + private static PartialFactorResult128 PartialFactorWide(in UInt128 value, in PrimeOrderCalculatorConfig config) + { + if (value <= UInt128.One) + { + return PartialFactorResult128.Empty; + } + + Dictionary counts = new(capacity: 8); + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + UInt128 remaining; + + // if (IsGpuPow2Allowed && value <= ulong.MaxValue) + // { + // Dictionary narrowCounts = new(capacity: 8); + // if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining)) + // { + // foreach (KeyValuePair entry in narrowCounts) + // { + // counts[(UInt128)entry.Key] = entry.Value; + // } + + // remaining = (UInt128)narrowRemaining; + // } + // else + // { + // remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); + // } + // } + // else + // { + remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); + // } + + List pending = new(); + if (remaining > UInt128.One) + { + pending.Add(remaining); + } + + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + Stack stack = new(); + stack.Push(remaining); + pending.Clear(); + + long timestamp = 0L; // reused for deadline checks. + while (stack.Count > 0) + { + UInt128 composite = stack.Pop(); + if (composite <= UInt128.One) + { + continue; + } + + if (IsPrimeWide(composite)) + { + AddFactor(counts, composite, 1); + continue; + } + + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + pending.Add(composite); + continue; + } + + if (!TryPollardRhoWide(composite, deadlineTimestamp, out UInt128 factor)) + { + pending.Add(composite); + continue; + } + + UInt128 quotient = composite / factor; + stack.Push(factor); + stack.Push(quotient); + } + } + + UInt128 cofactor = UInt128.One; + int pendingCount = pending.Count; + for (int i = 0; i < pendingCount; i++) + { + UInt128 composite = pending[i]; + if (IsPrimeWide(composite)) + { + AddFactor(counts, composite, 1); + } + else + { + BigInteger product = (BigInteger)cofactor * (BigInteger)composite; + cofactor = (UInt128)product; + } + } + + if (counts.Count == 0 && cofactor == value) + { + return new PartialFactorResult128(null, value, false, 0); + } + + ArrayPool? pool = null; + FactorEntry128[]? rented = counts.Count > 0 ? (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null; + int index = 0; + if (rented is not null) + { + foreach (KeyValuePair entry in counts) + { + rented[index] = new FactorEntry128(entry.Key, entry.Value); + index++; + } + } + + FactorEntry128[]? resultArray = null; + if (pool is not null) + { + Span span = rented.AsSpan(0, index); + span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + resultArray = new FactorEntry128[index]; + span.CopyTo(resultArray); + pool.Return(rented!, clearArray: false); + } + + bool fullyFactored = cofactor == UInt128.One; + return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index); + } + + private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary 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(in UInt128 n, long deadlineTimestamp, out UInt128 factor) + { + factor = UInt128.Zero; + if ((n & UInt128.One) == UInt128.Zero) + { + factor = (UInt128)2UL; + return true; + } + + UInt128 c = UInt128.One; + UInt128 x = (UInt128)2UL; + UInt128 y = x; + long timestamp = 0L; // reused for deadline checks. + + while (true) + { + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + return false; + } + + x = AdvancePolynomialWide(x, c, n); + y = AdvancePolynomialWide(y, c, n); + y = AdvancePolynomialWide(y, c, n); + + UInt128 diff = x > y ? x - y : y - x; + UInt128 d = BinaryGcdWide(diff, n); + if (d == n) + { + c += UInt128.One; + x = (UInt128)2UL; + y = x; + continue; + } + + if (d > UInt128.One) + { + factor = d; + return true; + } + } + } + + private static UInt128 AdvancePolynomialWide(in UInt128 x, in UInt128 c, in UInt128 modulus) + { + BigInteger value = (BigInteger)x; + value = (value * value + (BigInteger)c) % (BigInteger)modulus; + return (UInt128)value; + } + + private static UInt128 BinaryGcdWide(UInt128 a, UInt128 b) + { + if (a == UInt128.Zero) + { + return b; + } + + if (b == UInt128.Zero) + { + return a; + } + + int shift = TrailingZeroCountWide(a | b); + a >>= TrailingZeroCountWide(a); + + while (true) + { + b >>= TrailingZeroCountWide(b); + if (a > b) + { + (a, b) = (b, a); + } + + b -= a; + if (b == UInt128.Zero) + { + return a << shift; + } + } + } + + private static int TrailingZeroCountWide(in UInt128 value) + { + if (value == UInt128.Zero) + { + return 128; + } + + ulong low = (ulong)value; + if (low != 0UL) + { + return BitOperations.TrailingZeroCount(low); + } + + ulong high = (ulong)(value >> 64); + return 64 + BitOperations.TrailingZeroCount(high); + } + + private static void FactorCompletelyWide(in UInt128 value, Dictionary counts) + { + if (value <= UInt128.One) + { + return; + } + + if (IsPrimeWide(value)) + { + AddFactor(counts, value, 1); + return; + } + + UInt128 factor = PollardRhoWide(value); + FactorCompletelyWide(factor, counts); + FactorCompletelyWide(value / factor, counts); + } + + private static UInt128 PollardRhoWide(in UInt128 n) + { + if ((n & UInt128.One) == UInt128.Zero) + { + return (UInt128)2UL; + } + + UInt128 c = UInt128.One; + while (true) + { + UInt128 x = (UInt128)2UL; + UInt128 y = x; + UInt128 d = UInt128.One; + + while (d == UInt128.One) + { + x = AdvancePolynomialWide(x, c, n); + y = AdvancePolynomialWide(y, c, n); + y = AdvancePolynomialWide(y, c, n); + + UInt128 diff = x > y ? x - y : y - x; + d = BinaryGcdWide(diff, n); + } + + if (d != n) + { + return d; + } + + c += UInt128.One; + } + } + + private static bool IsPrimeWide(in UInt128 value) + { + if (value <= ulong.MaxValue) + { + // return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None); + return PrimeTester.IsPrimeCpu((ulong)value); + } + + if ((value & UInt128.One) == UInt128.Zero) + { + return false; + } + + uint[] smallPrimes = PrimesGenerator.SmallPrimes; + int length = smallPrimes.Length; + for (int i = 0; i < length; i++) + { + uint prime = smallPrimes[i]; + UInt128 primeValue = prime; + UInt128 primeSquare = primeValue * primeValue; + if (primeSquare > value) + { + break; + } + + if ((value % primeValue) == UInt128.Zero) + { + return false; + } + } + + BigInteger n = (BigInteger)value; + BigInteger d = n - BigInteger.One; + int s = 0; + while ((d & 1) == 0) + { + d >>= 1; + s++; + } + + int[] bases = new[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 }; + for (int i = 0; i < bases.Length; i++) + { + int baseValue = bases[i]; + if (baseValue <= 1) + { + continue; + } + + BigInteger a = baseValue; + if (a >= n) + { + continue; + } + + BigInteger x = BigInteger.ModPow(a, d, n); + if (x == BigInteger.One || x == n - BigInteger.One) + { + continue; + } + + bool witnessFound = true; + for (int r = 1; r < s; r++) + { + x = BigInteger.ModPow(x, 2, n); + if (x == n - BigInteger.One) + { + witnessFound = false; + break; + } + } + + if (witnessFound) + { + return false; + } + } + + return true; + } + + private static void AddFactor(Dictionary counts, in UInt128 prime, int exponent) + { + if (counts.TryGetValue(prime, out int existing)) + { + counts[prime] = existing + exponent; + } + else + { + counts[prime] = exponent; + } + } + + private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary) + { + public int Primary { get; } = primary; + + public bool SecondaryDescending { get; } = secondaryDescending; + + public UInt128 Secondary { get; } = secondary; + + public bool TertiaryDescending { get; } = tertiaryDescending; + + public UInt128 Tertiary { get; } = tertiary; + } + + private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count) + { + public FactorEntry128[]? Factors { get; } = factors; + + public UInt128 Cofactor { get; } = cofactor; + + public bool FullyFactored { get; } = fullyFactored; + + public int Count { get; } = count; + + public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0); + + public PartialFactorResult128 WithAdditionalPrime(UInt128 prime) + { + if (Factors is null) + { + FactorEntry128[] local = new FactorEntry128[1]; + local[0] = new FactorEntry128(prime, 1); + return new PartialFactorResult128(local, UInt128.One, true, 1); + } + + FactorEntry128[] extended = new FactorEntry128[Count + 1]; + Array.Copy(Factors, extended, Count); + extended[Count] = new FactorEntry128(prime, 1); + return new PartialFactorResult128(extended, UInt128.One, true, Count + 1); + } + } + + public static ulong Calculate( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong prime, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config) + { + if (prime <= 3UL) + { + return prime == 3UL ? 2UL : 1UL; + } + + ulong phi = prime - 1UL; + + #if DEVICE_GPU + if (PrimeOrderGpuHeuristics.TryCalculateOrder(gpu, prime, previousOrder, config, divisorData, out ulong gpuOrder)) + { + return gpuOrder; + } + + PartialFactorResult phiFactors = PartialFactor(gpu, phi, divisorData, config); + #elif DEVICE_HYBRID + PartialFactorResult phiFactors = PartialFactor(gpu, phi, divisorData, config); + #else + PartialFactorResult phiFactors = PartialFactor(phi, divisorData, config); + #endif + + ulong result; + if (phiFactors.Factors is null) + { + result = CalculateByFactorization(prime, divisorData); + phiFactors.Dispose(); + return result; + } + + #if DEVICE_HYBRID || DEVICE_GPU + result = RunHeuristicPipeline(gpu, prime, previousOrder, config, divisorData, phi, phiFactors); + #else + result = RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors); + #endif + phiFactors.Dispose(); + return result; + } + + public static UInt128 Calculate( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128? previousOrder, + in PrimeOrderCalculatorConfig config) + { + if (prime <= ulong.MaxValue) + { + ulong? previous = null; + if (previousOrder.HasValue) + { + UInt128 previousValue = previousOrder.Value; + previous = previousValue <= ulong.MaxValue ? (ulong)previousValue : ulong.MaxValue; + } + + ulong prime64 = (ulong)prime; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime64); + #if DEVICE_HYBRID || DEVICE_GPU + ulong order64 = Calculate(gpu, prime64, previous, divisorData, config); + #else + ulong order64 = Calculate(prime64, previous, divisorData, config); + #endif + return order64 == 0UL ? UInt128.Zero : (UInt128)order64; + } + + #if DEVICE_GPU + return CalculateWideInternal(gpu, prime, previousOrder, MontgomeryDivisorData.Empty, config); + #else + return CalculateWideInternal(prime, previousOrder, MontgomeryDivisorData.Empty, config); + #endif + } + + public static UInt128 Calculate( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config) + { + if (prime <= ulong.MaxValue) + { + ulong? previous = null; + if (previousOrder.HasValue) + { + UInt128 previousValue = previousOrder.Value; + previous = previousValue <= ulong.MaxValue ? (ulong)previousValue : ulong.MaxValue; + } + + ulong prime64 = (ulong)prime; + MontgomeryDivisorData effectiveDivisorData = divisorData.Equals(MontgomeryDivisorData.Empty) ? MontgomeryDivisorData.FromModulus(prime64) : divisorData; + #if DEVICE_HYBRID || DEVICE_GPU + ulong order64 = Calculate(gpu, prime64, previous, effectiveDivisorData, config); + #else + ulong order64 = Calculate(prime64, previous, effectiveDivisorData, config); + #endif + return order64 == 0UL ? UInt128.Zero : (UInt128)order64; + } + + #if DEVICE_GPU + return CalculateWideInternal(gpu, prime, previousOrder, divisorData, config); + #else + return CalculateWideInternal(prime, previousOrder, divisorData, config); + #endif + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static ulong CalculateByFactorization(ulong prime, in MontgomeryDivisorData divisorData) + { + ulong phi = prime - 1UL; + Dictionary counts = ThreadStaticPools.RentFactorCountDictionary(); + FactorCompletely(phi, counts, false); + PrintStats(counts.Count); + if (counts.Count == 0) + { + ThreadStaticPools.ReturnFactorCountDictionary(counts); + return phi; + } + + KeyValuePair[] entries = [.. counts]; + Array.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); + + ulong order = phi; + int entryCount = entries.Length; + + for (int i = 0; i < entryCount; i++) + { + ulong primeFactor = entries[i].Key; + int exponent = entries[i].Value; + for (int iteration = 0; iteration < exponent; iteration++) + { + ulong remainder = order.ReduceCycleRemainder(primeFactor); + if (remainder != 0UL) + { + break; + } + + // TODO: Benchmark performance with the division vs reusing remainder + ulong candidate = order - primeFactor * remainder; + // ulong candidate = order / primeFactor; + + if (candidate.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL) + { + order = candidate; + continue; + } + + break; + } + } + + ThreadStaticPools.ReturnFactorCountDictionary(counts); + return order; + } + + private static UInt128 CalculateWideInternal( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config) + { + if (prime <= UInt128.One) + { + return UInt128.One; + } + + if (prime == UInt128Numbers.Three) + { + return UInt128Numbers.Two; + } + + UInt128 phi = prime - UInt128.One; + PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); + if (phiFactors.Factors is null) + { + #if DEVICE_GPU + return FinishStrictlyWide(gpu, prime, divisorData); + #else + return FinishStrictlyWide(prime, divisorData); + #endif + } + + #if DEVICE_GPU + UInt128 candidateOrder = InitializeStartingOrderWide(gpu, prime, phi, divisorData); + candidateOrder = ExponentLoweringWide(gpu, candidateOrder, prime, divisorData, phiFactors); + #else + UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData); + candidateOrder = ExponentLoweringWide(candidateOrder, prime, divisorData, phiFactors); + #endif + + #if DEVICE_GPU + if (TryConfirmOrderWide(gpu, prime, candidateOrder, divisorData, config)) + #else + if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config)) + #endif + { + return candidateOrder; + } + + if (config.StrictMode) + { + #if DEVICE_GPU + return FinishStrictlyWide(gpu, prime, divisorData); + #else + return FinishStrictlyWide(prime, divisorData); + #endif + } + + #if DEVICE_GPU + if (TryHeuristicFinishWide(gpu, prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) + #else + if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) + #endif + { + return order; + } + + return candidateOrder; + } + + private static bool CheckCandidateViolation( + Span buffer, + ulong primeFactor, + int exponent, + ulong candidate, + ulong prime, + ref int powUsed, + int powBudget, + ref ExponentRemainderStepperCpu stepper) + { + ulong working = candidate; + int actual = 0; + while (actual < exponent) + { + if (working < primeFactor) + { + break; + } + + ulong reduced = working / primeFactor; + if (reduced == 0UL) + { + break; + } + + buffer[actual] = reduced; + working = reduced; + actual++; + } + + if (actual == 0) + { + return false; + } + + Span candidates = buffer[..actual]; + + stepper.Reset(); + int last = actual - 1; + + bool enforceBudget = powBudget > 0; + if (!enforceBudget) + { + powUsed++; + if (stepper.InitializeCpuIsUnity(candidates[last])) + { + return true; + } + + for (int j = last - 1; j >= 0; j--) + { + powUsed++; + if (stepper.ComputeNextIsUnity(candidates[j])) + { + return true; + } + } + + return false; + } + + if (powUsed >= powBudget) + { + return true; + } + + powUsed++; + if (stepper.InitializeCpuIsUnity(candidates[last])) + { + return true; + } + + for (int j = last - 1; j >= 0; j--) + { + if (powUsed >= powBudget) + { + return true; + } + + powUsed++; + if (stepper.ComputeNextIsUnity(candidates[j])) + { + return true; + } + } + + return false; + } + +#if DEVICE_CPU || DEVICE_HYBRID + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void FactorCompletely(ulong value, Dictionary counts, bool knownComposite) + { + if (value <= 1UL) + { + return; + } + + bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(value); + if (!knownComposite && isPrime) + { + AddFactor(counts, value, 1); + return; + } + + ulong factor = PollardRhoStrict(value); + ulong quotient = value / factor; + + isPrime = PrimeTesterByLastDigit.IsPrimeCpu(factor); + if (isPrime) + { + int exponent = 1; + ulong remaining = quotient; + while (remaining.ReduceCycleRemainder(factor) == 0UL) + { + remaining /= factor; + exponent++; + } + + AddFactor(counts, factor, exponent); + + if (remaining > 1UL) + { + FactorCompletely(remaining, counts, knownComposite: false); + } + + return; + } + + FactorCompletely(factor, counts, knownComposite: true); + + if (quotient == factor) + { + isPrime = false; + } + else + { + isPrime = PrimeTesterByLastDigit.IsPrimeCpu(quotient); + } + + if (isPrime) + { + AddFactor(counts, quotient, 1); + } + else + { + FactorCompletely(quotient, counts, knownComposite: true); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static void AddFactor(Dictionary counts, ulong prime, int exponent) + { + if (counts.TryGetValue(prime, out int existing)) + { + counts[prime] = existing + exponent; + } + else + { + counts[prime] = exponent; + } + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void AddFactorToCollector( + Span primes, + Span exponents, + ref int count, + ulong prime) + { + for (int i = 0; i < count; i++) + { + if (primes[i] == prime) + { + exponents[i]++; + return; + } + } + + primes[count] = prime; + exponents[count] = 1; + count++; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static void AddFactorToDictionary(Dictionary counts, ulong prime) + { + if (counts.TryGetValue(prime, out int existing)) + { + counts[prime] = existing + 1; + } + else + { + counts[prime] = 1; + } + } +#endif + +#if DEVICE_CPU || DEVICE_HYBRID + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void PopulateSmallPrimeFactors( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) + { + factorCount = 0; + if (value <= 1UL) + { + remaining = value; + return; + } + + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + + for (int i = 0; i < primeCount && value > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > limit || squares[i] > value) + { + break; + } + + ulong primeValue = primeCandidate; + ulong quotient = value / primeValue; + if (value - (quotient * primeValue) != 0UL) + { + continue; + } + + int exponent = ExtractSmallPrimeExponent(ref value, primeValue, quotient); + + primeTargets[factorCount] = primeValue; + exponentTargets[factorCount] = exponent; + factorCount++; + } + + remaining = value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static int ExtractSmallPrimeExponent(ref ulong value, ulong primeValue, ulong quotient) + { + ulong dividend = quotient; + int exponent = 1; + + if (dividend < primeValue) + { + value = dividend; + return exponent; + } + + quotient = dividend / primeValue; + ulong product = quotient * primeValue; + if (dividend - product != 0UL) + { + value = dividend; + return exponent; + } + + dividend = quotient; + exponent++; + + if (dividend < primeValue) + { + value = dividend; + return exponent; + } + + quotient = dividend / primeValue; + product = quotient * primeValue; + if (dividend - product != 0UL) + { + value = dividend; + return exponent; + } + + dividend = quotient; + exponent++; + + while (true) + { + quotient = dividend / primeValue; + product = quotient * primeValue; + + if (dividend - product != 0UL) + { + break; + } + + dividend = quotient; + exponent++; + } + + value = dividend; + return exponent; + } +#endif + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static void AddFactorToDictionary(Dictionary counts, UInt128 prime) + { + if (counts.TryGetValue(prime, out int existing)) + { + counts[prime] = existing + 1; + } + else + { + counts[prime] = 1; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) + { + return unchecked((x * x) % modulus + c) % modulus; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static ulong AdvancePolynomialTwice(ulong x, ulong c, ulong modulus) + { + x = AdvancePolynomial(x, c, modulus); + return AdvancePolynomial(x, c, modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void CollectFactors( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + Span primeSlots, + Span exponentSlots, + ref int factorCount, + List pending, + FixedCapacityStack compositeStack, + long deadlineTimestamp, + out bool pollardRhoDeadlineReached) + { + bool limitReached = false; + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + #if DEVICE_GPU + bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + #elif DEVICE_HYBRID + bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + #else + bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(composite); + #endif + + if (isPrime) + { + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); + continue; + } + + limitReached = Stopwatch.GetTimestamp() > deadlineTimestamp; + if (limitReached) + { + pending.Add(new (composite, knownComposite: true)); + continue; + } + + #if DEVICE_GPU + if (!TryPollardRho(gpu, composite, out ulong factor)) + #else + if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) + #endif + { + pending.Add(new (composite, knownComposite: true)); + continue; + } + + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + + pollardRhoDeadlineReached = limitReached; + } + +#if DEVICE_GPU || DEVICE_HYBRID + private static void CollectFactors( + PrimeOrderCalculatorAccelerator gpu, + ref Span primeSlots, + ref Span exponentSlots, + ref ulong[] factors, + ref int[] exponents, + ref int factorCount, + List pending, + FixedCapacityStack compositeStack, + long deadlineTimestamp, + ref bool pollardRhoDeadlineReached) + { + CollectFactors(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out bool reached); + pollardRhoDeadlineReached |= reached; + } +#endif + + +#if DEVICE_HYBRID || DEVICE_GPU + private static bool EvaluateSpecialMaxCandidates( + PrimeOrderCalculatorAccelerator gpu, + ReadOnlySpan factors, + ulong phi, + ulong prime, + in MontgomeryDivisorData divisorData) + { + int factorCount = factors.Length; + gpu.EnsureUlongInputOutputCapacity(factorCount); + + int acceleratorIndex = gpu.AcceleratorIndex; + ArrayView specialMaxFactorsView = gpu.InputView; + ArrayView specialMaxResultView = gpu.OutputUlongView2; + var kernelLauncher = gpu.SpecialMaxKernelLauncher; + + var stream = AcceleratorStreamPool.Rent(acceleratorIndex); + specialMaxFactorsView.SubView(0, factorCount).CopyFromCPU(stream, factors); + + kernelLauncher( + stream, + 1, + phi, + specialMaxFactorsView, + factorCount, + divisorData.Modulus, + specialMaxResultView); + + ulong result = 0UL; + specialMaxResultView.CopyToCPU(stream, ref result, 1); + stream.Synchronize(); + + AcceleratorStreamPool.Return(acceleratorIndex, stream); + return result != 0; + } +#endif + +#if DEVICE_CPU || DEVICE_HYBRID + private static bool EvaluateSpecialMaxCandidates( + Span buffer, + ReadOnlySpan factors, + ulong phi, + ulong prime, + in MontgomeryDivisorData divisorData) + { + int factorCount = factors.Length; + for (int i = 0; i < factorCount; i++) + { + buffer[i] = phi / factors[i]; + } + + buffer.Sort(); + + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + + if (stepper.InitializeCpuIsUnity(buffer[0])) + { + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + return false; + } + + for (int i = 1; i < factorCount; i++) + { + if (stepper.ComputeNextIsUnity(buffer[i])) + { + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + return false; + } + } + + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + return true; + } +#endif + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static UInt128 ExponentLoweringWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + UInt128 order, + in UInt128 prime, + in MontgomeryDivisorData divisorData, + in PartialFactorResult128 factors) + { + ArrayPool pool = ThreadStaticPools.FactorEntry128Pool; + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + // TODO: This should never trigger from production code - check + // if (length == 0) + // { + // return order; + // } + FactorEntry128[] tempArray = pool.Rent(length + 1); + Span buffer = tempArray.AsSpan(0, length); + factorSpan.CopyTo(buffer); + if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor)) + { + buffer[length] = new FactorEntry128(factors.Cofactor, 1); + length++; + } + + buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value)); + + for (int i = 0; i < length; i++) + { + UInt128 primeFactor = buffer[i].Value; + int exponent = buffer[i].Exponent; + for (int iteration = 0; iteration < exponent; iteration++) + { + #if DEVICE_GPU + if ((order % primeFactor) == UInt128.Zero) + #else + if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero) + #endif + { + UInt128 reduced = order / primeFactor; + #if DEVICE_GPU + if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One) + #else + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) + #endif + { + order = reduced; + continue; + } + } + + break; + } + } + + pool.Return(tempArray, clearArray: false); + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static UInt128 FinishStrictlyWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in MontgomeryDivisorData divisorData) + { + UInt128 phi = prime - UInt128.One; + FixedCapacityStack> uInt128IntDictionaryPool = ThreadStaticPools.UInt128IntDictionaryPool; + Dictionary counts = uInt128IntDictionaryPool.Rent(capacity: 8); + FactorCompletelyWide(phi, counts); + int entryCount = counts.Count; + PrintStats(entryCount); + if (entryCount == 0) + { + uInt128IntDictionaryPool.Return(counts); + return phi; + } + + FixedCapacityStack>> keyValuePairUInt128IntegerPool = ThreadStaticPools.KeyValuePairUInt128IntegerPool; + var entries = keyValuePairUInt128IntegerPool.Rent(entryCount); + entries.AddRange(counts); + uInt128IntDictionaryPool.Return(counts); + + entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); + + UInt128 order = phi; + for (int i = 0; i < entryCount; i++) + { + UInt128 primeFactor = entries[i].Key; + int exponent = entries[i].Value; + for (int iteration = 0; iteration < exponent; iteration++) + { + #if DEVICE_GPU + if ((order % primeFactor) != UInt128.Zero) + #else + if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero) + #endif + { + UInt128 candidate = order / primeFactor; + #if DEVICE_GPU + if (Pow2ModWide(gpu, candidate, prime, divisorData) == UInt128.One) + #else + if (Pow2ModWide(candidate, prime, divisorData) == UInt128.One) + #endif + { + order = candidate; + } + else + { + break; + } + } + } + } + + keyValuePairUInt128IntegerPool.Return(entries); + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + #if DEVICE_HYBRID || DEVICE_GPU + private static ulong InitializeStartingOrder(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong phi, in MontgomeryDivisorData divisorData) + #else + private static ulong InitializeStartingOrder(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) + #endif + { + ulong order = phi; + if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) + { + ulong half = phi >> 1; + #if DEVICE_HYBRID + if (Pow2EqualsOne(gpu, half, prime, divisorData)) + #else + if (Pow2EqualsOne(half, divisorData)) + #endif + { + order = half; + } + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static UInt128 InitializeStartingOrderWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128 phi, + in MontgomeryDivisorData divisorData) + { + UInt128 order = phi; + UInt128 mod8 = prime & (UInt128)7UL; + if (mod8 == UInt128.One || mod8 == (UInt128)7UL) + { + UInt128 half = phi >> 1; + #if DEVICE_GPU + if (Pow2ModWide(gpu, half, prime, divisorData) == UInt128.One) + #else + if (Pow2ModWide(half, prime, divisorData) == UInt128.One) + #endif + { + order = half; + } + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static bool Pow2EqualsOne(ulong exponent, in MontgomeryDivisorData divisorData) + { + return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; + } + +#if DEVICE_HYBRID || DEVICE_GPU + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static bool Pow2EqualsOne( + PrimeOrderCalculatorAccelerator gpu, + ulong exponent, + in MontgomeryDivisorData divisorData) + { + // TODO: swap to GPU pow2 when available; fallback to CPU for now. + return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static bool Pow2EqualsOne( + PrimeOrderCalculatorAccelerator gpu, + ulong exponent, + ulong prime, + in MontgomeryDivisorData divisorData) + { + return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; + } +#endif + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static ulong ExponentLowering( + ulong order, + PartialFactorResult factors, + in MontgomeryDivisorData divisorData) + { + int factorCount = factors.Count; + int newFactorCount = factorCount + 1; + Span factorsSpan = new(factors.Factors, 0, newFactorCount); + Span exponentsSpan = new(factors.Exponents, 0, newFactorCount); + + if (!factors.FullyFactored && factors.Cofactor > 1UL && factors.CofactorIsPrime) + { + factorsSpan[factorCount] = factors.Cofactor; + exponentsSpan[factorCount] = 1; + factors.Count = factorCount = newFactorCount; + factors.Cofactor = 1UL; + factors.FullyFactored = true; + factors.CofactorIsPrime = false; + } + + SortFactorExponentSpans(factorsSpan, exponentsSpan, factorCount); + + ExponentRemainderStepperCpu stepper = factors.ExponentRemainderStepper; + + Span stackCandidates = new(factors.StackCandidates); + Span stackEvaluations = new(factors.StackEvaluations); + + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factorsSpan[i]; + int exponent = exponentsSpan[i]; + ProcessExponentLoweringPrime(stackCandidates[..exponent], stackEvaluations[..exponent], ref order, primeFactor, exponent, ref stepper); + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void SortFactorExponentSpans(Span factorsSpan, Span exponentsSpan, int length) + { + if (length <= 1) + { + return; + } + + for (int i = 1; i < length; i++) + { + int current = i; + int previous = current - 1; + while (current > 0 && factorsSpan[current] < factorsSpan[previous]) + { + (factorsSpan[current], factorsSpan[previous]) = (factorsSpan[previous], factorsSpan[current]); + (exponentsSpan[current], exponentsSpan[previous]) = (exponentsSpan[previous], exponentsSpan[current]); + current--; + previous--; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void ProcessExponentLoweringPrime( + Span candidateBuffer, + Span evaluationBuffer, + ref ulong order, + ulong primeFactor, + int exponent, + ref ExponentRemainderStepperCpu stepper) + { + ulong working = order; + int actual = 0; + while (actual < exponent) + { + if (working < primeFactor) + { + break; + } + + ulong reduced = working / primeFactor; + if (reduced == 0UL) + { + break; + } + + candidateBuffer[actual] = reduced; + working = reduced; + actual++; + } + + if (actual == 0) + { + return; + } + + Span candidates = candidateBuffer[..actual]; + Span evaluations = evaluationBuffer[..actual]; + + stepper.Reset(); + + int last = actual - 1; + evaluations[last] = stepper.InitializeCpuIsUnity(candidates[last]); + for (int j = last - 1; j >= 0; j--) + { + evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); + } + + for (int j = 0; j < actual; j++) + { + if (!evaluations[j]) + { + break; + } + + order = candidates[j]; + } + } + +#if DEVICE_GPU + private static PartialFactorResult PartialFactor( + PrimeOrderCalculatorAccelerator gpu, + ulong value, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + // stackalloc is faster than pooling + + PartialFactorResult result = PartialFactorResult.Rent(divisorData); + ulong[] factors = result.Factors; + Span primeSlots = new(factors); + int[] exponents = result.Exponents; + Span exponentSlots = new(exponents); + + // TODO: Check if this is true for GPU path + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); + + List pending = result.PendingFactors; + pending.Clear(); + + Span exponentBuffer = exponentSlots[..PrimeOrderConstants.ExponentSlotLimit]; + Span primeBuffer = primeSlots[..PrimeOrderConstants.PrimeSlotLimit]; + Span compositeSlots = stackalloc ulong[PrimeOrderConstants.MaxCompositeSlots]; + Span factorScratch = stackalloc ulong[PrimeOrderConstants.ScratchSlots]; + + FixedCapacityStack compositeStack = result.CompositeStack; + compositeStack.Clear(); + + ulong working = value; + bool pollardRhoDeadlineReached = false; + int factorCount; + CollectSmallFactors(gpu, config, exponentBuffer, primeBuffer, ref working, ref pollardRhoDeadlineReached, out factorCount); + + Span primeSlotsBuffer = primeSlots[..factorCount]; + Span exponentSlotsBuffer = exponentSlots[..factorCount]; + primeBuffer[..factorCount].CopyTo(primeSlotsBuffer); + exponentBuffer[..factorCount].CopyTo(exponentSlotsBuffer); + + if (working > 1UL && !pollardRhoDeadlineReached) + { + compositeStack.Push(working); + long deadlineTimestamp = Stopwatch.GetTimestamp() + (config.MaxSeconds * Stopwatch.Frequency); + CollectFactors(gpu, ref primeSlotsBuffer, ref exponentSlotsBuffer, ref factors, ref exponents, ref factorCount, pending, compositeStack, deadlineTimestamp, ref pollardRhoDeadlineReached); + } + + primeSlotsBuffer = primeSlots[..factorCount]; + exponentSlotsBuffer = exponentSlots[..factorCount]; + primeSlotsBuffer.Sort(exponentSlotsBuffer); + + result.Count = factorCount; + exponentSlotsBuffer.CopyTo(result.Exponents); + result.PendingFactors.Clear(); + result.PendingFactors.AddRange(pending); + result.Cofactor = pollardRhoDeadlineReached ? working : 1UL; + result.CofactorIsPrime = false; + + return result; + } +#else + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static PartialFactorResult PartialFactor( + #if DEVICE_HYBRID + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong value, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + // stackalloc is faster than pooling + + PartialFactorResult result = PartialFactorResult.Rent(divisorData); + ulong[] factors = result.Factors; + Span primeSlots = new(factors); + int[] exponents = result.Exponents; + Span exponentSlots = new(exponents); + + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); + + List pending = result.PendingFactors; + pending.Clear(); + + PopulateSmallPrimeFactors( + value, + config.SmallFactorLimit, + primeSlots, + exponentSlots, + out int factorCount, + out value); + + bool cofactorContainsComposite; + if (value > 1UL) + { + if (config.PollardRhoMilliseconds > 0) + { + FixedCapacityStack compositeStack = result.CompositeStack; + compositeStack.Clear(); + compositeStack.Push(value); + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + + #if DEVICE_HYBRID + CollectFactors(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out cofactorContainsComposite); + #else + CollectFactors(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out cofactorContainsComposite); + #endif + + if (cofactorContainsComposite) + { + while (compositeStack.Count > 0) + { + pending.Add(new(compositeStack.Pop(), knownComposite: false)); + } + } + } + else + { + pending.Add(new (value, knownComposite: false)); + } + } + + ulong cofactor = 1UL; + bool isPrime; + int pendingCount = pending.Count; + cofactorContainsComposite = false; + + for (int index = 0; index < pendingCount; index++) + { + PartialFactorPendingEntry entry = pending[index]; + // Reuse "value" parameter for testing if a number is composite to limit registry pressure. + value = entry.Value; + + if (entry.KnownComposite) + { + cofactor = checked(cofactor * value); + cofactorContainsComposite = true; + continue; + } + + if (!entry.HasKnownPrimality) + { + #if DEVICE_HYBRID + isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(value) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, value); + #else + isPrime = PrimeTesterByLastDigit.IsPrimeCpu(value); + #endif + + entry.WithPrimality(isPrime); + pending[index] = entry; + } + + if (entry.IsPrime) + { + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, value); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * value); + cofactorContainsComposite = true; + } + } + + if (cofactorContainsComposite || cofactor <= 1UL) + { + isPrime = false; + } + else + { + #if DEVICE_HYBRID + isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(cofactor) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor); + #else + isPrime = PrimeTesterByLastDigit.IsPrimeCpu(cofactor); + #endif + } + + // This will never happen in production code. We'll always get at least 1 factor + // if (factorCount == 0) + // { + // if (cofactor == value) + // (...) + // } + + Array.Sort(factors, exponents, 0, factorCount); + + result.Cofactor = cofactor; + result.FullyFactored = cofactor == 1UL; + result.Count = factorCount; + result.CofactorIsPrime = isPrime; + + return result; + } +#endif + + private static UInt128 Pow2ModWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 exponent, + in UInt128 modulus, + in MontgomeryDivisorData divisorData) + { + if (modulus == UInt128.One) + { + return UInt128.Zero; + } + + #if DEVICE_GPU + if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue) + { + ulong prime64 = (ulong)modulus; + ulong exponent64 = (ulong)exponent; + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent64, prime64, out ulong remainder, divisorData); + if (status == GpuPow2ModStatus.Success) + { + return remainder; + } + } + + GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, modulus, out UInt128 wideRemainder); + if (wideStatus == GpuPow2ModStatus.Success) + { + return wideRemainder; + } + + #endif + return exponent.Pow2MontgomeryModWindowed(modulus); + } + + private static ulong RunHeuristicPipeline( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong prime, + ulong? previousOrder, + in PrimeOrderCalculatorConfig config, + in MontgomeryDivisorData divisorData, + ulong phi, + PartialFactorResult phiFactors) + { + #if DEVICE_HYBRID || DEVICE_GPU + if (phiFactors.FullyFactored && TrySpecialMax(gpu, phi, prime, phiFactors, divisorData)) + #else + if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData)) + #endif + { + return phi; + } + + #if DEVICE_HYBRID || DEVICE_GPU + ulong candidateOrder = InitializeStartingOrder(gpu, prime, phi, divisorData); + #else + ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData); + #endif + candidateOrder = ExponentLowering(candidateOrder, phiFactors, divisorData); + + #if DEVICE_HYBRID || DEVICE_GPU + if (TryConfirmOrder(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) + #else + if (TryConfirmOrder(prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) + #endif + { + return candidateOrder; + } + + if (config.StrictMode) + { + orderFactors?.Dispose(); + return CalculateByFactorization(prime, divisorData); + } + + #if DEVICE_HYBRID || DEVICE_GPU + if (TryHeuristicFinish(gpu, prime, candidateOrder, previousOrder, divisorData, config, orderFactors, out ulong order)) + #else + if (TryHeuristicFinish(prime, candidateOrder, previousOrder, divisorData, config, orderFactors, out ulong order)) + #endif + { + return order; + } + + return candidateOrder; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static UInt128 RunHeuristicPipelineWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + in UInt128 phi, + in PartialFactorResult128 phiFactors) + { + #if DEVICE_GPU + if (phiFactors.FullyFactored && TrySpecialMaxWide(gpu, phi, prime, divisorData, phiFactors)) + #else + if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors)) + #endif + { + return phi; + } + + #if DEVICE_GPU + UInt128 candidateOrder = InitializeStartingOrderWide(gpu, prime, phi, divisorData); + candidateOrder = ExponentLoweringWide(gpu, candidateOrder, prime, divisorData, phiFactors); + + if (TryConfirmOrderWide(gpu, prime, candidateOrder, divisorData, config)) + #else + UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData); + candidateOrder = ExponentLoweringWide(candidateOrder, prime, divisorData, phiFactors); + + if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config)) + #endif + { + return candidateOrder; + } + + if (config.StrictMode) + { + #if DEVICE_GPU + return FinishStrictlyWide(gpu, prime, divisorData); + #else + return FinishStrictlyWide(prime, divisorData); + #endif + } + + #if DEVICE_GPU + if (TryHeuristicFinishWide(gpu, prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) + #else + if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) + #endif + { + return order; + } + + return candidateOrder; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) + { + ulong previous = previousOrder ?? 0UL; + bool hasPrevious = previousOrder.HasValue; + ulong threshold1 = prime >> 3; + ulong threshold2 = prime >> 2; + ulong threshold3 = (prime * 3UL) >> 3; + int previousGroup = hasPrevious ? GetGroup(previous, threshold1, threshold2, threshold3) : 1; + + candidates.Sort((x, y) => CompareCandidates(x, y, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3)); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static int CompareCandidates(ulong left, ulong right, ulong previous, int previousGroup, bool hasPrevious, ulong threshold1, ulong threshold2, ulong threshold3) + { + int leftGroup = GetGroup(left, threshold1, threshold2, threshold3); + int rightGroup = GetGroup(right, threshold1, threshold2, threshold3); + if (leftGroup == 0) + { + return rightGroup == 0 ? 0 : 1; + } + + if (rightGroup == 0) + { + return -1; + } + + bool leftIsGe = !hasPrevious || left >= previous; + bool rightIsGe = !hasPrevious || right >= previous; + int rightPrimary = ComputePrimary(rightGroup, rightIsGe, previousGroup); + + int compare = ComputePrimary(leftGroup, leftIsGe, previousGroup); + compare = compare.CompareTo(rightPrimary); + if (compare != 0) + { + return compare; + } + + ulong leftDistance = hasPrevious ? (left > previous ? left - previous : previous - left) : left; + ulong rightDistance = hasPrevious ? (right > previous ? right - previous : previous - right) : right; + + bool leftGroupDescending = leftGroup == 3; + if (leftDistance == rightDistance) + { + return leftGroupDescending ? (left > right ? -1 : (left < right ? 1 : 0)) : left.CompareTo(right); + } + + if (leftGroupDescending) + { + return leftDistance > rightDistance ? -1 : 1; + } + + return leftDistance.CompareTo(rightDistance); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static int GetGroup(ulong value, ulong threshold1, ulong threshold2, ulong threshold3) + { + if (value <= threshold1) + { + return 1; + } + + if (value <= threshold2) + { + return 2; + } + + if (value <= threshold3) + { + return 3; + } + + return 0; + } + + private static void SortCandidatesWide(UInt128 prime, in UInt128? previousOrder, List candidates) + { + UInt128 previous = previousOrder ?? UInt128.Zero; + int previousGroup = previousOrder.HasValue ? GetGroupWide(previousOrder.Value, prime) : 1; + bool hasPrevious = previousOrder.HasValue; + + candidates.Sort((x, y) => + { + CandidateKey128 keyX = BuildKeyWide(x, prime, previous, previousGroup, hasPrevious); + CandidateKey128 keyY = BuildKeyWide(y, prime, previous, previousGroup, hasPrevious); + + int primary = keyX.Primary.CompareTo(keyY.Primary); + if (primary != 0) + { + return primary; + } + + int secondary = CompareComponents(keyX.SecondaryDescending, keyX.Secondary, keyY.SecondaryDescending, keyY.Secondary); + if (secondary != 0) + { + return secondary; + } + + return CompareComponents(keyX.TertiaryDescending, keyX.Tertiary, keyY.TertiaryDescending, keyY.Tertiary); + }); + } + + private static bool TryConfirmCandidate( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong prime, + ulong candidate, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + ref int powUsed, + int powBudget) + { + #if DEVICE_HYBRID || DEVICE_GPU + PartialFactorResult factorization = PartialFactor(gpu, candidate, divisorData, config); + #else + PartialFactorResult factorization = PartialFactor(candidate, divisorData, config); + #endif + + if (factorization.Factors is null) + { + factorization.Dispose(); + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) + { + factorization.Dispose(); + return false; + } + + factorization.WithAdditionalPrime(factorization.Cofactor); + } + + ReadOnlySpan factorsSpan = factorization.Factors; + ReadOnlySpan exponentsSpan = factorization.Exponents; + int length = factorization.Count; + + const int StackExponentCapacity = 32; + const int ExponentHardLimit = 256; + + Span stackBuffer = stackalloc ulong[StackExponentCapacity]; + ulong[]? heapCandidateArray = null; + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + bool violates = false; + + for (int i = 0; i < length; i++) + { + ulong primeFactor = factorsSpan[i]; + int exponent = exponentsSpan[i]; + + if (exponent > ExponentHardLimit) + { + throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); + } + + if (exponent <= StackExponentCapacity) + { + if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) + { + violates = true; + break; + } + + continue; + } + + heapCandidateArray ??= FixedCapacityPools.ExclusiveUlongArray.Rent(ExponentHardLimit); + if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) + { + violates = true; + break; + } + } + + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + if (heapCandidateArray is not null) + { + FixedCapacityPools.ExclusiveUlongArray.Return(heapCandidateArray); + } + + factorization.Dispose(); + return !violates; + } + + private static bool TryConfirmCandidateWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128 candidate, + in PrimeOrderCalculatorConfig config, + ref int powUsed, + int powBudget, + in MontgomeryDivisorData divisorData) + { + PartialFactorResult128 factorization = PartialFactorWide(candidate, config); + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= UInt128.One) + { + return false; + } + + if (!IsPrimeWide(factorization.Cofactor)) + { + return false; + } + + factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + } + + ReadOnlySpan span = factorization.Factors; + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + UInt128 primeFactor = span[i].Value; + UInt128 reduced = candidate; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + #if DEVICE_GPU + if ((reduced % primeFactor) != UInt128.Zero) + #else + if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) + #endif + { + break; + } + + reduced /= primeFactor; + if (powUsed >= powBudget && powBudget > 0) + { + return false; + } + + powUsed++; + #if DEVICE_GPU + if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One) + #else + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) + #endif + { + return false; + } + } + } + + return true; + } + + private static bool TryConfirmOrder( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong prime, + ulong candidateOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + out PartialFactorResult? orderFactors) + { + orderFactors = null; + ArgumentOutOfRangeException.ThrowIfEqual(candidateOrder, 0UL); + + #if DEVICE_HYBRID + if (!Pow2EqualsOne(gpu, candidateOrder, prime, divisorData)) + #else + if (!Pow2EqualsOne(candidateOrder, divisorData)) + #endif + { + return false; + } + + #if DEVICE_HYBRID || DEVICE_GPU + PartialFactorResult factorization = PartialFactor(gpu, candidateOrder, divisorData, config); + #else + PartialFactorResult factorization = PartialFactor(candidateOrder, divisorData, config); + #endif + if (factorization.Factors is null) + { + factorization.Dispose(); + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) + { + factorization.Dispose(); + return false; + } + + factorization.WithAdditionalPrime(factorization.Cofactor); + } + + if (!ValidateOrderAgainstFactors(prime, candidateOrder, divisorData, factorization)) + { + orderFactors = factorization; + return false; + } + + factorization.Dispose(); + return true; + } + + private static bool TryConfirmOrderWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128 order, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config) + { + if (order == UInt128.Zero) + { + return false; + } + + #if DEVICE_GPU + if (Pow2ModWide(gpu, order, prime, divisorData) != UInt128.One) + #else + if (Pow2ModWide(order, prime, divisorData) != UInt128.One) + #endif + { + return false; + } + + PartialFactorResult128 factorization = PartialFactorWide(order, config); + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= UInt128.One) + { + return false; + } + + if (!IsPrimeWide(factorization.Cofactor)) + { + return false; + } + + factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + } + + ReadOnlySpan span = factorization.Factors; + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + UInt128 primeFactor = span[i].Value; + UInt128 reduced = order; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) + { + break; + } + + reduced /= primeFactor; + #if DEVICE_GPU + if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One) + #else + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) + #endif + { + return false; + } + } + } + + return true; + } + + private static bool TryHeuristicFinish( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong prime, + ulong candidateOrder, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + PartialFactorResult? orderFactors, + out ulong order) + { + order = 0UL; + if (candidateOrder <= 1UL) + { + orderFactors?.Dispose(); + return false; + } + + #if DEVICE_GPU || DEVICE_HYBRID + PartialFactorResult localOrderFactors = orderFactors ?? PartialFactor(gpu, candidateOrder, divisorData, config); + #else + PartialFactorResult localOrderFactors = orderFactors ?? PartialFactor(candidateOrder, divisorData, config); + #endif + try + { + if (localOrderFactors.Factors is null) + { + return false; + } + + if (!localOrderFactors.FullyFactored) + { + if (localOrderFactors.Cofactor <= 1UL || !localOrderFactors.CofactorIsPrime) + { + return false; + } + + localOrderFactors.WithAdditionalPrime(localOrderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; + List candidates = ThreadStaticPools.RentUlongList(capacity); + candidates.Clear(); + ulong[] factorArray = localOrderFactors.Factors!; + int[] exponentsArray = localOrderFactors.Exponents!; + + BuildCandidates(candidateOrder, factorArray, exponentsArray, localOrderFactors.Count, candidates, capacity); + int candidateCount = candidates.Count; + if (candidateCount == 0) + { + ThreadStaticPools.ReturnUlongList(candidates); + return false; + } + + SortCandidates(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; + int powUsed = 0; + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + bool powStepperInitialized = false; + + int index = 0; + #if DEVICE_GPU || DEVICE_HYBRID + bool allowGpuBatch = true; + Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; + ArrayPool pool = ThreadStaticPools.UlongPool; + #endif + + while (index < candidateCount && powUsed < powBudget) + { + int remaining = candidateCount - index; + int budgetRemaining = powBudget - powUsed; + int batchSize = budgetRemaining <= 0 ? remaining : Math.Min(remaining, budgetRemaining); + + #if DEVICE_GPU || DEVICE_HYBRID + batchSize = Math.Min(batchSize, PrimeOrderConstants.MaxGpuBatchSize); + #endif + + if (batchSize <= 0) + { + break; + } + + ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); + + #if DEVICE_GPU || DEVICE_HYBRID + ulong[]? gpuPool = null; + Span pooledGpuRemainders = default; + bool gpuSuccess = false; + bool gpuStackRemainders = false; + GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; + + if (allowGpuBatch) + { + if (batchSize <= PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount) + { + Span localRemainders = stackGpuRemainders[..batchSize]; + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, localRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + gpuSuccess = true; + gpuStackRemainders = true; + } + } + else + { + gpuPool = pool.Rent(batchSize); + Span pooledRemainders = gpuPool.AsSpan(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, pooledRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + pooledGpuRemainders = pooledRemainders; + gpuSuccess = true; + } + else + { + pool.Return(gpuPool, clearArray: false); + gpuPool = null; + } + } + + if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) + { + allowGpuBatch = false; + } + } + #endif + + for (int i = 0; i < batchSize && powUsed < powBudget; i++) + { + ulong candidate = batch[i]; + + #if DEVICE_GPU || DEVICE_HYBRID + bool equalsOne; + if (gpuSuccess) + { + powUsed++; + ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; + equalsOne = remainderValue == 1UL; + } + else + #endif + { + powUsed++; + if (!powStepperInitialized) + { + powStepperInitialized = true; + if (!powStepper.InitializeCpuIsUnity(candidate)) + { + continue; + } + } + else if (!powStepper.ComputeNextIsUnity(candidate)) + { + continue; + } + #if DEVICE_GPU || DEVICE_HYBRID + equalsOne = true; + #endif + } + + if (powUsed > powBudget && powBudget > 0) + { + break; + } + + #if DEVICE_GPU || DEVICE_HYBRID + if (!TryConfirmCandidate(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) + #else + if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget)) + #endif + { + continue; + } + + #if DEVICE_GPU || DEVICE_HYBRID + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + #endif + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + order = candidate; + return true; + } + + #if DEVICE_GPU || DEVICE_HYBRID + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + #endif + + index += batchSize; + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + + return false; + } + finally + { + localOrderFactors.Dispose(); + } + } + + private static bool TryHeuristicFinishWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 prime, + in UInt128 order, + in UInt128? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + out UInt128 result) + { + result = UInt128.Zero; + if (order <= UInt128.One) + { + return false; + } + + PartialFactorResult128 orderFactors = PartialFactorWide(order, config); + if (orderFactors.Factors is null) + { + return false; + } + + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= UInt128.One) + { + return false; + } + + if (!IsPrimeWide(orderFactors.Cofactor)) + { + return false; + } + + orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; + FixedCapacityStack> uInt128ListPool = ThreadStaticPools.UInt128ListPool; + List candidates = uInt128ListPool.Rent(capacity); + FactorEntry128[] factorArray = orderFactors.Factors!; + + BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); + + int candidateCount = candidates.Count; + if (candidateCount == 0) + { + uInt128ListPool.Return(candidates); + return false; + } + + SortCandidatesWide(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; + int powUsed = 0; + + for (int i = 0; i < candidateCount; i++) + { + if (powUsed >= powBudget && powBudget > 0) + { + break; + } + + UInt128 candidate = candidates[i]; + powUsed++; + + #if DEVICE_GPU + if (Pow2ModWide(gpu, candidate, prime, divisorData) != UInt128.One) + #else + if (Pow2ModWide(candidate, prime, divisorData) != UInt128.One) + #endif + { + continue; + } + + #if DEVICE_GPU + if (!TryConfirmCandidateWide(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData)) + #else + if (!TryConfirmCandidateWide(prime, candidate, config, ref powUsed, powBudget, divisorData)) + #endif + { + continue; + } + + result = candidate; + uInt128ListPool.Return(candidates); + return true; + } + + uInt128ListPool.Return(candidates); + return false; + } + +#if DEVICE_GPU + private static bool TryPollardRho( + PrimeOrderCalculatorAccelerator gpu, + ulong n, + out ulong factor) + { + Span randomStateSpan = stackalloc ulong[1]; + randomStateSpan[0] = ThreadStaticDeterministicRandomGpu.Exclusive.State; + + Span factoredSpan = stackalloc byte[1]; + Span factorSpan = stackalloc ulong[1]; + + int acceleratorIndex = gpu.AcceleratorIndex; + var kernel = gpu.PollardRhoKernel; + + AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); + gpu.InputView.CopyFromCPU(stream, randomStateSpan); + + kernel.Launch(stream, 1, n, 1, gpu.InputView, gpu.OutputByteView, gpu.OutputUlongView); + + gpu.OutputByteView.CopyToCPU(stream, factoredSpan); + gpu.OutputUlongView.CopyToCPU(stream, factorSpan); + gpu.InputView.CopyToCPU(stream, randomStateSpan); + stream.Synchronize(); + + AcceleratorStreamPool.Return(acceleratorIndex, stream); + ThreadStaticDeterministicRandomGpu.Exclusive.SetState(randomStateSpan[0]); + + bool factored = factoredSpan[0] != 0; + factor = factored ? factorSpan[0] : 0UL; + return factored; + } +#else + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TryPollardRho( + ulong n, + long deadlineTimestamp, + out ulong factor) + { + factor = 0UL; + if ((n & 1UL) == 0UL) + { + factor = 2UL; + return true; + } + + long timestamp; // reused for deadline checks. + while (true) + { + ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + return false; + } + + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomialTwice(y, c, n); + ulong diff = x > y ? x - y : y - x; + d = diff.BinaryGcd(n); + } + + if (d != n) + { + factor = d; + return true; + } + } + } +#endif + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static long CreateDeadlineTimestamp(int milliseconds) + { + if (milliseconds <= 0) + { + return long.MaxValue; + } + + long ticksPerMillisecond = Stopwatch.Frequency / 1000; + return Stopwatch.GetTimestamp() + (ticksPerMillisecond * milliseconds); + } + +#if DEVICE_HYBRID + private static int _cpuCount; + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool RunOnCpu() + { + int gpuRatio = EnvironmentConfiguration.GpuRatio; + if (gpuRatio == 0) + { + return false; + } + + if (gpuRatio == 1) + { + while (true) + { + if (Interlocked.CompareExchange(ref _cpuCount, 0, 1) == 0) + { + return true; + } + + if (Interlocked.CompareExchange(ref _cpuCount, 1, 0) == 1) + { + return false; + } + + Thread.Yield(); + } + } + + int cpuCount = Interlocked.Increment(ref _cpuCount); + if (cpuCount == gpuRatio) + { + Interlocked.Add(ref _cpuCount, -gpuRatio); + } + else if (cpuCount > gpuRatio) + { + cpuCount -= gpuRatio; + } + + return cpuCount != gpuRatio; + } +#endif + + private static bool TrySpecialMax( + #if DEVICE_HYBRID || DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + ulong phi, + ulong prime, + PartialFactorResult factors, + in MontgomeryDivisorData divisorData) + { + int length = factors.Count; + ReadOnlySpan factorSpan = new(factors.Factors, 0, length); + Span stackBuffer = stackalloc ulong[length]; + + #if DEVICE_HYBRID + if (length <= 8) + { + return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData); + } + + return EvaluateSpecialMaxCandidates(gpu, factorSpan, phi, prime, divisorData); + #elif DEVICE_GPU + return EvaluateSpecialMaxCandidates(gpu, factorSpan, phi, prime, divisorData); + #else + return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData); + #endif + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TrySpecialMaxWide( + #if DEVICE_GPU + PrimeOrderCalculatorAccelerator gpu, + #endif + in UInt128 phi, + in UInt128 prime, + in MontgomeryDivisorData divisorData, + in PartialFactorResult128 factors) + { + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + for (int i = 0; i < length; i++) + { + UInt128 factor = factorSpan[i].Value; + UInt128 reduced = phi / factor; + #if DEVICE_GPU + if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One) + #else + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) + #endif + { + return false; + } + } + + return true; + } + +#if DEVICE_GPU + private const int GpuSmallPrimeFactorSlots = 64; + + // GPU specific. It doesn't make sense to consolidate into one method. + private static bool TryPopulateSmallPrimeFactors( + PrimeOrderCalculatorAccelerator gpu, + ulong value, + uint limit, + Dictionary counts, + out int factorCount, + out ulong remaining) + { + var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); + var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); + Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + remaining = value; + + gpu.EnsureSmallPrimeFactorSlotsCapacity(GpuSmallPrimeFactorSlots); + int acceleratorIndex = gpu.AcceleratorIndex; + + var kernelLauncher = gpu.SmallPrimeFactorKernelLauncher; + ArrayView1D smallPrimeFactorCountSlotView = gpu.OutputIntView2; + ArrayView1D smallPrimeFactorRemainingSlotView = gpu.InputView; + ArrayView1D smallPrimeFactorPrimeSlotsView = gpu.OutputUlongView; + ArrayView1D smallPrimeFactorExponentSlotsView = gpu.OutputIntView; + ArrayView1D smallPrimeFactorsSquaresView = gpu.SmallPrimeFactorSquares; + ArrayView1D smallPrimeFactorsPrimesView = gpu.SmallPrimeFactorPrimes; + + var stream = AcceleratorStreamPool.Rent(acceleratorIndex); + kernelLauncher( + stream, + 1, + value, + limit, + smallPrimeFactorsPrimesView, + smallPrimeFactorsSquaresView, + (int)smallPrimeFactorsPrimesView.Length, + smallPrimeFactorPrimeSlotsView, + smallPrimeFactorExponentSlotsView, + smallPrimeFactorCountSlotView, + smallPrimeFactorRemainingSlotView); + + factorCount = 0; + smallPrimeFactorCountSlotView.CopyToCPU(stream, ref factorCount, 1); + stream.Synchronize(); + + factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); + + if (factorCount > 0) + { + primeBuffer = primeBuffer[..factorCount]; + smallPrimeFactorPrimeSlotsView.CopyToCPU(stream, primeBuffer); + + exponentBuffer = exponentBuffer[..factorCount]; + smallPrimeFactorExponentSlotsView.CopyToCPU(stream, exponentBuffer); + } + + smallPrimeFactorRemainingSlotView.CopyToCPU(stream, ref remaining, 1); + stream.Synchronize(); + + AcceleratorStreamPool.Return(acceleratorIndex, stream); + + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + counts.Add(primeValue, exponent); + } + + ThreadStaticPools.UlongPool.Return(primeBufferArray, clearArray: false); + ThreadStaticPools.IntPool.Return(exponentBufferArray, clearArray: false); + return true; + } + + // GPU specific. It doesn't make sense to consolidate into one method. + private static void CollectSmallFactors( + PrimeOrderCalculatorAccelerator gpu, + in PrimeOrderCalculatorConfig config, + Span exponentBuffer, + Span primeBuffer, + ref ulong working, + ref bool pollardRhoDeadlineReached, + out int factorCount) + { + Dictionary counts = ThreadStaticPools.UlongIntDictionaryPool.Rent(capacity: 64); + counts.Clear(); + + bool populated = TryPopulateSmallPrimeFactors(gpu, working, (uint)PrimeNumberConstants.GpuSmallPrimeBound, counts, out factorCount, out ulong remaining); + pollardRhoDeadlineReached = false; + if (!populated) + { + factorCount = 0; + working = remaining; + ThreadStaticPools.UlongIntDictionaryPool.Return(counts); + return; + } + + if (remaining > 1UL) + { + if (counts.Count < GpuSmallPrimeFactorSlots) + { + counts.Add(remaining, 1); + working = 1UL; + } + else + { + pollardRhoDeadlineReached = true; + working = remaining; + } + } + + int collectedCount = 0; + foreach (KeyValuePair pair in counts) + { + if (collectedCount >= PrimeOrderConstants.ExponentSlotLimit) + { + pollardRhoDeadlineReached = true; + break; + } + + primeBuffer[collectedCount] = pair.Key; + exponentBuffer[collectedCount] = pair.Value; + collectedCount++; + } + + factorCount = collectedCount; + ThreadStaticPools.UlongIntDictionaryPool.Return(counts); + } + +#endif + + private static bool ValidateOrderAgainstFactors( + ulong prime, + ulong order, + in MontgomeryDivisorData divisorData, + PartialFactorResult factorization) + { + ReadOnlySpan factorsSpan = factorization.Factors!; + ReadOnlySpan exponentsSpan = factorization.Exponents!; + + int length = factorization.Count; + + const int StackExponentCapacity = 32; + + Span stackBuffer = stackalloc ulong[StackExponentCapacity]; + + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + bool violates = false; + + for (int i = 0; i < length; i++) + { + ulong primeFactor = factorsSpan[i]; + int exponent = exponentsSpan[i]; + + if (ValidateOrderForFactor(stackBuffer[..exponent], primeFactor, exponent, order, ref stepper)) + { + violates = true; + break; + } + } + + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + + return !violates; + } + + private static bool ValidateOrderForFactor( + Span buffer, + ulong primeFactor, + int exponent, + ulong order, + ref ExponentRemainderStepperCpu stepper) + { + ulong working = order; + int actual = 0; + while (actual < exponent) + { + if (working < primeFactor) + { + break; + } + + ulong reduced = working / primeFactor; + if (reduced == 0UL) + { + break; + } + + buffer[actual] = reduced; + working = reduced; + actual++; + } + + if (actual == 0) + { + return false; + } + + Span candidates = buffer[..actual]; + + stepper.Reset(); + + int last = actual - 1; + if (stepper.InitializeCpuIsUnity(candidates[last])) + { + return true; + } + + for (int j = last - 1; j >= 0; j--) + { + if (stepper.ComputeNextIsUnity(candidates[j])) + { + return true; + } + } + + return false; + } +} diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index c6af7d6..a3ccc2c 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -100,6 +100,31 @@ namespace PerfectNumbers.Core [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public static void Return(this FixedCapacityStack> pool, Dictionary dictionary) => pool.Push(dictionary); + [ThreadStatic] + private static FixedCapacityStack>? _ulongIntDictionaryPool; + + public static FixedCapacityStack> UlongIntDictionaryPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + get => _ulongIntDictionaryPool ??= new(PerfectNumberConstants.DefaultPoolCapacity); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static Dictionary Rent(this FixedCapacityStack> pool, int capacity) + { + if (pool.Pop() is { } pooled) + { + pooled.EnsureCapacity(capacity); + pooled.Clear(); + return pooled; + } + + return new(capacity); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static void Return(this FixedCapacityStack> pool, Dictionary dictionary) => pool.Push(dictionary); + [ThreadStatic] private static ArrayPool? _factorEntryPool; diff --git a/PerfectNumbers.Core/UInt128Numbers.cs b/PerfectNumbers.Core/UInt128Numbers.cs index 3ed08be..f94a331 100644 --- a/PerfectNumbers.Core/UInt128Numbers.cs +++ b/PerfectNumbers.Core/UInt128Numbers.cs @@ -5,6 +5,7 @@ namespace PerfectNumbers.Core public static class UInt128Numbers { public static readonly UInt128 Two = 2UL; + public static readonly UInt128 Three = 3UL; public static readonly UInt128 Four = 4UL; public static readonly UInt128 SixtyFour = 64UL; public static readonly UInt128 OneHundredTwentySeven = 127UL; diff --git a/PerfectNumbers.Core/UnboundedTaskScheduler.cs b/PerfectNumbers.Core/UnboundedTaskScheduler.cs index 33fcfbf..fee3bf2 100644 --- a/PerfectNumbers.Core/UnboundedTaskScheduler.cs +++ b/PerfectNumbers.Core/UnboundedTaskScheduler.cs @@ -28,7 +28,7 @@ public sealed class UnboundedTaskScheduler : TaskScheduler public static void ConfigureThreadCount(int threadCount) { _configuredThreadCount = threadCount; - PrimeOrderCalculator.PrimeOrderCalculatorConfig.ConfigureHeuristicDefault(threadCount); + PrimeOrderCalculatorConfig.ConfigureHeuristicDefault(threadCount); // On EvenPerfectBitScanner's execution path the scheduler instance is created after configuration, // so updating an existing pool would never run. Leave the update path disabled for now.