From d3695a3dc3ca09fd8a15cca7ffc777d217d2325b Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 15 Oct 2025 02:17:13 +0200 Subject: [PATCH] Pools made thread static. Various small optimizations. "in" parameter added to read-only method parameters --- ...rsenneNumberDivisorByDivisorTesterTests.cs | 4 +- ...MersenneNumberDivisorByDivisorCpuTester.cs | 113 +- PerfectNumbers.Core/FactorEntry128.cs | 8 + ...MersenneNumberDivisorByDivisorGpuTester.cs | 4 +- .../IMersenneNumberDivisorByDivisorTester.cs | 4 +- .../MersenneNumberDivisorByDivisorTester.cs | 390 ++- PerfectNumbers.Core/PartialFactorResult.cs | 6 +- .../PrimeOrderCalculator.Wide.cs | 169 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 2749 ++++++++--------- PerfectNumbers.Core/ThreadStaticPools.cs | 41 + 10 files changed, 1713 insertions(+), 1775 deletions(-) create mode 100644 PerfectNumbers.Core/FactorEntry128.cs create mode 100644 PerfectNumbers.Core/ThreadStaticPools.cs diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs index 42b4068..9131db8 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs @@ -70,7 +70,7 @@ public class MersenneNumberDivisorByDivisorTesterTests return new DummySession(); } - public void PrepareCandidates(ReadOnlySpan primes, Span allowedMaxValues) + public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) { for (int i = 0; i < primes.Length; i++) { @@ -80,7 +80,7 @@ public class MersenneNumberDivisorByDivisorTesterTests private sealed class DummySession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession { - public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) + public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits) { throw new NotSupportedException(); } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 49baab4..e9fdd0c 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,18 +1,14 @@ -using System; using System.Collections.Concurrent; -using System.Collections.Generic; -using System.Numerics; -using PerfectNumbers.Core; +using System.Runtime.CompilerServices; namespace PerfectNumbers.Core.Cpu; public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester { private readonly object _sync = new(); - private readonly ConcurrentBag _sessionPool = new(); + private readonly ConcurrentBag _sessionPool = []; private ulong _divisorLimit; private ulong _lastStatusDivisor; - private bool _isConfigured; private int _batchSize = 1_024; @@ -24,55 +20,25 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv public void ConfigureFromMaxPrime(ulong maxPrime) { - lock (_sync) - { - _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - _lastStatusDivisor = 0UL; - _isConfigured = true; - } + _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); + _lastStatusDivisor = 0UL; } public ulong DivisorLimit - { + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] get { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return _divisorLimit; - } + return _divisorLimit; } } - public ulong GetAllowedMaxDivisor(ulong prime) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong GetAllowedMaxDivisor(ulong prime) => ComputeAllowedMaxDivisor(prime, _divisorLimit); + + public bool IsPrime(ulong prime, out bool divisorsExhausted) { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return ComputeAllowedMaxDivisor(prime, _divisorLimit); - } - } - - public bool IsPrime(ulong prime, out bool divisorsExhausted) - { - ulong allowedMax; - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - } + ulong allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); if (allowedMax < 3UL) { @@ -81,13 +47,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } ulong processedCount; - ulong lastProcessed; bool processedAll; bool composite = CheckDivisors( prime, allowedMax, - out lastProcessed, out processedAll, out processedCount); @@ -109,26 +73,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return true; } - public void PrepareCandidates(ReadOnlySpan primes, Span allowedMaxValues) + public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) { - if (allowedMaxValues.Length < primes.Length) - { - throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues)); - } - - ulong divisorLimit; - - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - divisorLimit = _divisorLimit; - } - - for (int index = 0; index < primes.Length; index++) + ulong divisorLimit = _divisorLimit; + int length = primes.Length; + for (int index = 0; index < length; index++) { allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit); } @@ -136,21 +85,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession() { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + if (_sessionPool.TryTake(out DivisorScanSession? session)) + { + session.Reset(); + return session; + } - if (_sessionPool.TryTake(out DivisorScanSession? session)) - { - session.Reset(); - return session; - } - - return new DivisorScanSession(this); - } + return new DivisorScanSession(this); } private void ReturnSession(DivisorScanSession session) @@ -161,11 +102,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool CheckDivisors( ulong prime, ulong allowedMax, - out ulong lastProcessed, out bool processedAll, out ulong processedCount) { - lastProcessed = 0UL; processedCount = 0UL; processedAll = false; @@ -214,7 +153,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv { ulong candidate = (ulong)divisor; processedCount++; - lastProcessed = candidate; bool admissible = lastIsSeven ? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9) @@ -348,13 +286,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _disposed = false; } - public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) + public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits) { - if (_disposed) - { - throw new ObjectDisposedException(nameof(DivisorScanSession)); - } - int length = primes.Length; if (length == 0) { diff --git a/PerfectNumbers.Core/FactorEntry128.cs b/PerfectNumbers.Core/FactorEntry128.cs new file mode 100644 index 0000000..ea7b978 --- /dev/null +++ b/PerfectNumbers.Core/FactorEntry128.cs @@ -0,0 +1,8 @@ +namespace PerfectNumbers.Core; + +public readonly struct FactorEntry128(in UInt128 value, int exponent) +{ + public readonly UInt128 Value = value; + + public readonly int Exponent = exponent; +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 2220711..4adde79 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -148,7 +148,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return true; } - public void PrepareCandidates(ReadOnlySpan primes, Span allowedMaxValues) + public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) { if (allowedMaxValues.Length < primes.Length) { @@ -583,7 +583,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _disposed = false; } - public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) + public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits) { if (_disposed) { diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index e700843..fee848f 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -16,10 +16,10 @@ public interface IMersenneNumberDivisorByDivisorTester IDivisorScanSession CreateDivisorSession(); - void PrepareCandidates(ReadOnlySpan primes, Span allowedMaxValues); + void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues); public interface IDivisorScanSession : IDisposable { - void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits); + void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits); } } diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index f5e07a9..efbeb16 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -10,129 +10,126 @@ namespace PerfectNumbers.Core; public static class MersenneNumberDivisorByDivisorTester { - public static void Run( - List candidates, - IMersenneNumberDivisorByDivisorTester tester, - Dictionary? previousResults, - ulong startPrime, - Action markComposite, - Action clearComposite, - Action printResult, - int threadCount) - { - if (candidates.Count == 0) - { - Console.WriteLine("No candidates were provided for --mersenne=bydivisor."); - return; - } + public static void Run( + List candidates, + IMersenneNumberDivisorByDivisorTester tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount) + { + if (candidates.Count == 0) + { + Console.WriteLine("No candidates were provided for --mersenne=bydivisor."); + return; + } - _ = threadCount; + _ = threadCount; - candidates.Sort(); - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - int candidateCount = candidateSpan.Length; + candidates.Sort(); + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + int candidateCount = candidateSpan.Length; - bool applyStartPrime = startPrime > 0UL; - int skippedByPreviousResults = 0; + bool applyStartPrime = startPrime > 0UL; + int skippedByPreviousResults = 0; - if (previousResults is not null && previousResults.Count > 0) - { - int recordedCount = previousResults.Count; - ulong[] recordedCandidates = ArrayPool.Shared.Rent(recordedCount); - try - { - Span recordedSpan = recordedCandidates.AsSpan(0, recordedCount); - int recordedIndex = 0; - foreach (ulong value in previousResults.Keys) - { - recordedSpan[recordedIndex++] = value; - } + ArrayPool pool = ThreadStaticPools.UlongPool; + if (previousResults is not null && previousResults.Count > 0) + { + int recordedCount = previousResults.Count; + ulong[] recordedCandidates = pool.Rent(recordedCount); - recordedSpan.Sort(); + Span recordedSpan = recordedCandidates.AsSpan(0, recordedCount); + int recordedIndex = 0; + foreach (ulong value in previousResults.Keys) + { + recordedSpan[recordedIndex++] = value; + } - int writeIndex = 0; - int skipIndex = 0; + recordedSpan.Sort(); - for (int readIndex = 0; readIndex < candidateCount; readIndex++) - { - ulong candidate = candidateSpan[readIndex]; + int writeIndex = 0; + int skipIndex = 0; - while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate) - { - skipIndex++; - } + for (int readIndex = 0; readIndex < candidateCount; readIndex++) + { + ulong candidate = candidateSpan[readIndex]; - if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate) - { - skippedByPreviousResults++; - skipIndex++; - continue; - } + while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate) + { + skipIndex++; + } - candidateSpan[writeIndex++] = candidate; - } + if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate) + { + skippedByPreviousResults++; + skipIndex++; + continue; + } - if (writeIndex < candidateCount) - { - candidates.RemoveRange(writeIndex, candidateCount - writeIndex); - candidateSpan = CollectionsMarshal.AsSpan(candidates); - candidateCount = candidateSpan.Length; - } - } - finally - { - ArrayPool.Shared.Return(recordedCandidates, clearArray: true); - } - } + candidateSpan[writeIndex++] = candidate; + } - if (skippedByPreviousResults > 0) - { - Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results."); - } + if (writeIndex < candidateCount) + { + candidates.RemoveRange(writeIndex, candidateCount - writeIndex); + candidateSpan = CollectionsMarshal.AsSpan(candidates); + candidateCount = candidateSpan.Length; + } - int startIndex = 0; - if (applyStartPrime) - { - startIndex = candidates.BinarySearch(startPrime); - if (startIndex < 0) - { - startIndex = ~startIndex; - } - else - { - while (startIndex > 0 && candidates[startIndex - 1] >= startPrime) - { - startIndex--; - } - } - } + pool.Return(recordedCandidates, clearArray: true); + } - List primesToTest = candidates; - Span primesSpan = CollectionsMarshal.AsSpan(primesToTest); - ulong maxPrime = 0UL; - int primeWriteIndex = 0; + if (skippedByPreviousResults > 0) + { + Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results."); + } + + int startIndex = 0; + if (applyStartPrime) + { + startIndex = candidates.BinarySearch(startPrime); + if (startIndex < 0) + { + startIndex = ~startIndex; + } + else + { + while (startIndex > 0 && candidates[startIndex - 1] >= startPrime) + { + startIndex--; + } + } + } + + List primesToTest = candidates; + Span primesSpan = CollectionsMarshal.AsSpan(primesToTest); + ulong maxPrime = 0UL; + int primeWriteIndex = 0; // This implementation is terribly slow, while this method expect prime p given as --filter-p input already. // We don't need to additionally check it. - // using IEnumerator primeEnumerator = Prime.Numbers.GetEnumerator(); - // bool hasPrime = primeEnumerator.MoveNext(); - // ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL; + // using IEnumerator primeEnumerator = Prime.Numbers.GetEnumerator(); + // bool hasPrime = primeEnumerator.MoveNext(); + // ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL; - for (int index = startIndex; index < candidateCount; index++) - { - ulong candidate = primesSpan[index]; + for (int index = startIndex; index < candidateCount; index++) + { + ulong candidate = primesSpan[index]; - if (candidate <= 1UL) - { - markComposite(); - printResult(candidate, false, false, false); - continue; - } + if (candidate <= 1UL) + { + markComposite(); + printResult(candidate, false, false, false); + continue; + } // This implementation is terribly slow, while this method expect prime p given as --filter-p input already. // We don't need to additionally check it. - + // while (hasPrime && currentPrime < candidate) // { // hasPrime = primeEnumerator.MoveNext(); @@ -154,123 +151,118 @@ public static class MersenneNumberDivisorByDivisorTester primesSpan[primeWriteIndex++] = candidate; - if (candidate > maxPrime) - { - maxPrime = candidate; - } - } + if (candidate > maxPrime) + { + maxPrime = candidate; + } + } - if (primeWriteIndex < primesToTest.Count) - { - primesToTest.RemoveRange(primeWriteIndex, primesToTest.Count - primeWriteIndex); - } + if (primeWriteIndex < primesToTest.Count) + { + primesToTest.RemoveRange(primeWriteIndex, primesToTest.Count - primeWriteIndex); + } - if (primesToTest.Count == 0) - { - if (applyStartPrime) - { - Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); - } - else - { - Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering."); - } + if (primesToTest.Count == 0) + { + if (applyStartPrime) + { + Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); + } + else + { + Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering."); + } - return; - } + return; + } - if (maxPrime <= 1UL) - { - Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor."); - return; - } + if (maxPrime <= 1UL) + { + Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor."); + return; + } - tester.ConfigureFromMaxPrime(maxPrime); + tester.ConfigureFromMaxPrime(maxPrime); - int primeCount = primesToTest.Count; - ulong[] primeBatch = ArrayPool.Shared.Rent(primeCount); - ulong[] allowedMaxBatch = ArrayPool.Shared.Rent(primeCount); - List filteredPrimes = new(primeCount); + int primeCount = primesToTest.Count; + ulong[] primeBatch = pool.Rent(primeCount); + ulong[] allowedMaxBatch = pool.Rent(primeCount); + List filteredPrimes = new(primeCount); - try - { - Span primeSpan = primeBatch.AsSpan(0, primeCount); - Span allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount); + Span primeSpan = primeBatch.AsSpan(0, primeCount); + Span allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount); - for (int i = 0; i < primeCount; i++) - { - primeSpan[i] = primesToTest[i]; - } + for (int i = 0; i < primeCount; i++) + { + primeSpan[i] = primesToTest[i]; + } - tester.PrepareCandidates(primeSpan, allowedMaxSpan); + tester.PrepareCandidates(primeSpan, allowedMaxSpan); - for (int i = 0; i < primeCount; i++) - { - ulong prime = primeSpan[i]; - ulong allowedMax = allowedMaxSpan[i]; + for (int i = 0; i < primeCount; i++) + { + ulong prime = primeSpan[i]; + ulong allowedMax = allowedMaxSpan[i]; - if (allowedMax < 3UL) - { - clearComposite(); - printResult(prime, true, true, true); - continue; - } + if (allowedMax < 3UL) + { + clearComposite(); + printResult(prime, true, true, true); + continue; + } - filteredPrimes.Add(prime); - } - } - finally - { - ArrayPool.Shared.Return(allowedMaxBatch, clearArray: true); - ArrayPool.Shared.Return(primeBatch, clearArray: true); - } + filteredPrimes.Add(prime); + } - if (filteredPrimes.Count == 0) - { - if (applyStartPrime) - { - Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); - } + pool.Return(allowedMaxBatch, clearArray: true); + pool.Return(primeBatch, clearArray: true); - return; - } + if (filteredPrimes.Count == 0) + { + if (applyStartPrime) + { + Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); + } - int workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount; - if (workerCount < 1) - { - workerCount = 1; - } + return; + } - void ProcessPrime(ulong prime) - { - bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted); + int workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount; + if (workerCount < 1) + { + workerCount = 1; + } - if (!isPrime) - { - markComposite(); - printResult(prime, true, true, false); - return; - } + void ProcessPrime(ulong prime) + { + bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted); - clearComposite(); - printResult(prime, true, divisorsExhausted, true); - } + if (!isPrime) + { + markComposite(); + printResult(prime, true, true, false); + return; + } - if (workerCount == 1) - { - foreach (ulong prime in filteredPrimes) - { - ProcessPrime(prime); - } - } - else - { - ParallelOptions options = new() - { - MaxDegreeOfParallelism = workerCount - }; + clearComposite(); + printResult(prime, true, divisorsExhausted, true); + } - Parallel.ForEach(filteredPrimes, options, ProcessPrime); - } - } + if (workerCount == 1) + { + foreach (ulong prime in filteredPrimes) + { + ProcessPrime(prime); + } + } + else + { + ParallelOptions options = new() + { + MaxDegreeOfParallelism = workerCount + }; + + Parallel.ForEach(filteredPrimes, options, ProcessPrime); + } + } } diff --git a/PerfectNumbers.Core/PartialFactorResult.cs b/PerfectNumbers.Core/PartialFactorResult.cs index dc9c1b0..c57665d 100644 --- a/PerfectNumbers.Core/PartialFactorResult.cs +++ b/PerfectNumbers.Core/PartialFactorResult.cs @@ -38,8 +38,8 @@ namespace PerfectNumbers.Core return instance; } - public static readonly PartialFactorResult Empty = new PartialFactorResult - { + public static readonly PartialFactorResult Empty = new() + { Cofactor = 1UL, FullyFactored = true, }; @@ -70,7 +70,7 @@ namespace PerfectNumbers.Core if (factors is not null) { Factors = null; - ArrayPool.Shared.Return(factors, clearArray: false); + ThreadStaticPools.FactorEntryPool.Return(factors, clearArray: false); } _next = s_poolHead; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs index fc06285..b1c89c9 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs @@ -10,10 +10,10 @@ namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { public static UInt128 Calculate( - UInt128 prime, - UInt128? previousOrder, - PrimeOrderSearchConfig config, - PrimeOrderHeuristicDevice device) + in UInt128 prime, + in UInt128? previousOrder, + in PrimeOrderSearchConfig config, + in PrimeOrderHeuristicDevice device) { var scope = UsePow2Mode(device); MontgomeryDivisorData divisorData; @@ -49,7 +49,7 @@ internal static partial class PrimeOrderCalculator return result; } - private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) + private static UInt128 CalculateWideInternal(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { if (prime <= UInt128.One) { @@ -79,12 +79,12 @@ internal static partial class PrimeOrderCalculator } private static UInt128 RunHeuristicPipelineWide( - UInt128 prime, - UInt128? previousOrder, + in UInt128 prime, + in UInt128? previousOrder, in MontgomeryDivisorData divisorData, - PrimeOrderSearchConfig config, - UInt128 phi, - PartialFactorResult128 phiFactors) + in PrimeOrderSearchConfig config, + in UInt128 phi, + in PartialFactorResult128 phiFactors) { if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors)) { @@ -112,12 +112,12 @@ internal static partial class PrimeOrderCalculator return candidateOrder; } - private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData) + private static UInt128 FinishStrictlyWide(in UInt128 prime, in MontgomeryDivisorData divisorData) { return CalculateByFactorizationWide(prime, divisorData); } - private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors) + private static bool TrySpecialMaxWide(in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors) { ReadOnlySpan factorSpan = factors.Factors; int length = factors.Count; @@ -134,7 +134,7 @@ internal static partial class PrimeOrderCalculator return true; } - private static UInt128 InitializeStartingOrderWide(UInt128 prime, UInt128 phi, in MontgomeryDivisorData divisorData) + private static UInt128 InitializeStartingOrderWide(in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData divisorData) { UInt128 order = phi; UInt128 mod8 = prime & (UInt128)7UL; @@ -150,19 +150,20 @@ internal static partial class PrimeOrderCalculator return order; } - private static UInt128 ExponentLoweringWide(UInt128 order, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors) + private static UInt128 ExponentLoweringWide(UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors) { - FactorEntry128[]? tempArray = null; + 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); try { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - if (length == 0) - { - return order; - } - tempArray = ArrayPool.Shared.Rent(length + 1); Span buffer = tempArray.AsSpan(0, length); factorSpan.CopyTo(buffer); @@ -198,14 +199,11 @@ internal static partial class PrimeOrderCalculator } finally { - if (tempArray is not null) - { - ArrayPool.Shared.Return(tempArray, clearArray: false); - } + pool.Return(tempArray, clearArray: false); } } - private static bool TryConfirmOrderWide(UInt128 prime, UInt128 order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) + private static bool TryConfirmOrderWide(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { if (order == UInt128.Zero) { @@ -263,11 +261,11 @@ internal static partial class PrimeOrderCalculator } private static bool TryHeuristicFinishWide( - UInt128 prime, - UInt128 order, - UInt128? previousOrder, + in UInt128 prime, + in UInt128 order, + in UInt128? previousOrder, in MontgomeryDivisorData divisorData, - PrimeOrderSearchConfig config, + in PrimeOrderSearchConfig config, out UInt128 result) { result = UInt128.Zero; @@ -339,7 +337,7 @@ internal static partial class PrimeOrderCalculator return false; } - private static void SortCandidatesWide(UInt128 prime, UInt128? previousOrder, List candidates) + 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; @@ -367,9 +365,9 @@ internal static partial class PrimeOrderCalculator } private static CandidateKey128 BuildKeyWide( - UInt128 value, - UInt128 prime, - UInt128 previous, + in UInt128 value, + in UInt128 prime, + in UInt128 previous, int previousGroup, bool hasPrevious) { @@ -391,7 +389,7 @@ internal static partial class PrimeOrderCalculator return new CandidateKey128(primary, false, distance, false, value); } - private static int CompareComponents(bool descendingX, UInt128 valueX, bool descendingY, UInt128 valueY) + private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY) { if (descendingX == descendingY) { @@ -416,7 +414,7 @@ internal static partial class PrimeOrderCalculator return valueX > valueY ? 1 : -1; } - private static int GetGroupWide(UInt128 value, UInt128 prime) + private static int GetGroupWide(in UInt128 value, in UInt128 prime) { UInt128 threshold1 = prime >> 3; if (value <= threshold1) @@ -439,7 +437,7 @@ internal static partial class PrimeOrderCalculator return 0; } - private static void BuildCandidatesWide(UInt128 order, FactorEntry128[] factors, int count, List candidates, int limit) + private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List candidates, int limit) { if (count == 0) { @@ -452,10 +450,10 @@ internal static partial class PrimeOrderCalculator } private static void BuildCandidatesRecursiveWide( - UInt128 order, - ReadOnlySpan factors, + in UInt128 order, + in ReadOnlySpan factors, int index, - UInt128 divisorProduct, + in UInt128 divisorProduct, List candidates, int limit) { @@ -511,7 +509,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmCandidateWide(UInt128 prime, UInt128 candidate, PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData) + private static bool TryConfirmCandidateWide(in UInt128 prime, in UInt128 candidate, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData) { PartialFactorResult128 factorization = PartialFactorWide(candidate, config); if (factorization.Factors is null) @@ -564,7 +562,7 @@ internal static partial class PrimeOrderCalculator return true; } - private static PartialFactorResult128 PartialFactorWide(UInt128 value, PrimeOrderSearchConfig config) + private static PartialFactorResult128 PartialFactorWide(in UInt128 value, in PrimeOrderSearchConfig config) { if (value <= UInt128.One) { @@ -664,7 +662,8 @@ internal static partial class PrimeOrderCalculator return new PartialFactorResult128(null, value, false, 0); } - FactorEntry128[]? rented = counts.Count > 0 ? ArrayPool.Shared.Rent(counts.Count) : null; + ArrayPool? pool = null; + FactorEntry128[]? rented = counts.Count > 0 ? (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null; int index = 0; if (rented is not null) { @@ -676,20 +675,20 @@ internal static partial class PrimeOrderCalculator } FactorEntry128[]? resultArray = null; - if (rented is not 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); - ArrayPool.Shared.Return(rented, clearArray: false); + pool.Return(rented!, clearArray: false); } bool fullyFactored = cofactor == UInt128.One; return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index); } - private static UInt128 PopulateSmallPrimeFactorsCpuWide(UInt128 value, uint limit, Dictionary counts) + private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary counts) { UInt128 remaining = value; uint[] primes = PrimesGenerator.SmallPrimes; @@ -729,7 +728,7 @@ internal static partial class PrimeOrderCalculator return remaining; } - private static bool TryPollardRhoWide(UInt128 n, Stopwatch stopwatch, long budgetTicks, out UInt128 factor) + private static bool TryPollardRhoWide(in UInt128 n, Stopwatch stopwatch, long budgetTicks, out UInt128 factor) { factor = UInt128.Zero; if ((n & UInt128.One) == UInt128.Zero) @@ -771,7 +770,7 @@ internal static partial class PrimeOrderCalculator } } - private static UInt128 AdvancePolynomialWide(UInt128 x, UInt128 c, UInt128 modulus) + private static UInt128 AdvancePolynomialWide(in UInt128 x, in UInt128 c, in UInt128 modulus) { BigInteger value = (BigInteger)x; value = (value * value + (BigInteger)c) % (BigInteger)modulus; @@ -809,7 +808,7 @@ internal static partial class PrimeOrderCalculator } } - private static int TrailingZeroCountWide(UInt128 value) + private static int TrailingZeroCountWide(in UInt128 value) { if (value == UInt128.Zero) { @@ -865,7 +864,7 @@ internal static partial class PrimeOrderCalculator return order; } - private static void FactorCompletelyWide(UInt128 value, Dictionary counts) + private static void FactorCompletelyWide(in UInt128 value, Dictionary counts) { if (value <= UInt128.One) { @@ -883,7 +882,7 @@ internal static partial class PrimeOrderCalculator FactorCompletelyWide(value / factor, counts); } - private static UInt128 PollardRhoWide(UInt128 n) + private static UInt128 PollardRhoWide(in UInt128 n) { if ((n & UInt128.One) == UInt128.Zero) { @@ -916,7 +915,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool IsPrimeWide(UInt128 value) + private static bool IsPrimeWide(in UInt128 value) { if (value <= ulong.MaxValue) { @@ -996,7 +995,7 @@ internal static partial class PrimeOrderCalculator return true; } - private static UInt128 Pow2ModWide(UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData) + private static UInt128 Pow2ModWide(in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData) { if (modulus == UInt128.One) { @@ -1026,7 +1025,7 @@ internal static partial class PrimeOrderCalculator return exponent.Pow2MontgomeryModWindowed(modulus); } - private static void AddFactor(Dictionary counts, UInt128 prime, int exponent) + private static void AddFactor(Dictionary counts, in UInt128 prime, int exponent) { if (counts.TryGetValue(prime, out int existing)) { @@ -1038,60 +1037,30 @@ internal static partial class PrimeOrderCalculator } } - private readonly struct CandidateKey128 - { - public CandidateKey128(int primary, bool secondaryDescending, UInt128 secondary, bool tertiaryDescending, UInt128 tertiary) - { - Primary = primary; - SecondaryDescending = secondaryDescending; - Secondary = secondary; - TertiaryDescending = tertiaryDescending; - Tertiary = tertiary; - } + private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary) + { + public int Primary { get; } = primary; - public int Primary { get; } + public bool SecondaryDescending { get; } = secondaryDescending; - public bool SecondaryDescending { get; } + public UInt128 Secondary { get; } = secondary; - public UInt128 Secondary { get; } + public bool TertiaryDescending { get; } = tertiaryDescending; - public bool TertiaryDescending { get; } + public UInt128 Tertiary { get; } = tertiary; + } - public UInt128 Tertiary { get; } - } + private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count) + { + public FactorEntry128[]? Factors { get; } = factors; - private readonly struct FactorEntry128 - { - public FactorEntry128(UInt128 value, int exponent) - { - Value = value; - Exponent = exponent; - } + public UInt128 Cofactor { get; } = cofactor; - public UInt128 Value { get; } + public bool FullyFactored { get; } = fullyFactored; - public int Exponent { get; } - } + public int Count { get; } = count; - private readonly struct PartialFactorResult128 - { - public PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count) - { - Factors = factors; - Cofactor = cofactor; - FullyFactored = fullyFactored; - Count = count; - } - - public FactorEntry128[]? Factors { get; } - - public UInt128 Cofactor { get; } - - public bool FullyFactored { get; } - - public int Count { get; } - - public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0); + public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0); public PartialFactorResult128 WithAdditionalPrime(UInt128 prime) { diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 1deac16..7008dad 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -10,1407 +10,1402 @@ namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - internal enum PrimeOrderMode - { - Heuristic, - Strict, - } + internal enum PrimeOrderMode + { + Heuristic, + Strict, + } - internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) - { - public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 10_240, maxPowChecks: 24, PrimeOrderMode.Heuristic); - public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); + internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) + { + public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 10_240, maxPowChecks: 24, PrimeOrderMode.Heuristic); + public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); - public readonly uint SmallFactorLimit = smallFactorLimit; - public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; - public readonly int MaxPowChecks = maxPowChecks; - public readonly PrimeOrderMode Mode = mode; - } + public readonly uint SmallFactorLimit = smallFactorLimit; + public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; + public readonly int MaxPowChecks = maxPowChecks; + public readonly PrimeOrderMode Mode = mode; + } - internal enum PrimeOrderHeuristicDevice - { - Cpu, - Gpu, - } + internal enum PrimeOrderHeuristicDevice + { + Cpu, + Gpu, + } - [ThreadStatic] - private static bool s_pow2ModeInitialized; + [ThreadStatic] + private static bool s_pow2ModeInitialized; - [ThreadStatic] - private static bool s_allowGpuPow2; + [ThreadStatic] + private static bool s_allowGpuPow2; - private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; + private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; - [ThreadStatic] - private static PrimeOrderHeuristicDevice s_deviceMode; + [ThreadStatic] + private static PrimeOrderHeuristicDevice s_deviceMode; - [ThreadStatic] - private static bool s_debugLoggingEnabled; + [ThreadStatic] + private static bool s_debugLoggingEnabled; - [ThreadStatic] - private static List? s_partialFactorPending; + [ThreadStatic] + private static List? s_partialFactorPending; - [ThreadStatic] - private static Stack? s_partialFactorCompositeStack; + [ThreadStatic] + private static Stack? s_partialFactorCompositeStack; - [ThreadStatic] - private static Dictionary? s_partialFactorCounts; + [ThreadStatic] + private static Dictionary? s_partialFactorCounts; - private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; + private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; - private readonly struct Pow2ModeScope - { - private readonly bool _previousInitialized; - private readonly bool _previousAllow; - private readonly PrimeOrderHeuristicDevice _previousDevice; + private readonly struct Pow2ModeScope + { + private readonly bool _previousInitialized; + private readonly bool _previousAllow; + private readonly PrimeOrderHeuristicDevice _previousDevice; - public Pow2ModeScope(PrimeOrderHeuristicDevice device) - { - _previousInitialized = s_pow2ModeInitialized; - _previousAllow = s_allowGpuPow2; - _previousDevice = s_deviceMode; - s_pow2ModeInitialized = true; - s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; - s_deviceMode = device; - } + public Pow2ModeScope(PrimeOrderHeuristicDevice device) + { + _previousInitialized = s_pow2ModeInitialized; + _previousAllow = s_allowGpuPow2; + _previousDevice = s_deviceMode; + s_pow2ModeInitialized = true; + s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; + s_deviceMode = device; + } - public void Dispose() - { - s_allowGpuPow2 = _previousAllow; - s_pow2ModeInitialized = _previousInitialized; - s_deviceMode = _previousDevice; - } - } + public void Dispose() + { + s_allowGpuPow2 = _previousAllow; + s_pow2ModeInitialized = _previousInitialized; + s_deviceMode = _previousDevice; + } + } - private readonly struct DebugLoggingScope - { - private readonly bool _previous; + private readonly struct DebugLoggingScope + { + private readonly bool _previous; - public DebugLoggingScope(bool enabled) - { - _previous = s_debugLoggingEnabled; - s_debugLoggingEnabled = enabled; - } + public DebugLoggingScope(bool enabled) + { + _previous = s_debugLoggingEnabled; + s_debugLoggingEnabled = enabled; + } - public void Dispose() - { - s_debugLoggingEnabled = _previous; - } - } + public void Dispose() + { + s_debugLoggingEnabled = _previous; + } + } - private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device); + private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device); - private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled); + private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled); - [Conditional("DEBUG")] - private static void DebugLog(string message) - { - if (s_debugLoggingEnabled) - { - Console.WriteLine(message); - } - } + [Conditional("DEBUG")] + private static void DebugLog(string message) + { + if (s_debugLoggingEnabled) + { + Console.WriteLine(message); + } + } - [Conditional("DEBUG")] - private static void DebugLog(Func messageFactory) - { - if (s_debugLoggingEnabled) - { - Console.WriteLine(messageFactory()); - } - } + [Conditional("DEBUG")] + private static void DebugLog(Func messageFactory) + { + if (s_debugLoggingEnabled) + { + Console.WriteLine(messageFactory()); + } + } - public static ulong Calculate( - ulong prime, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderSearchConfig config, - PrimeOrderHeuristicDevice device) - { - var scope = UsePow2Mode(device); - ulong order = CalculateInternal(prime, previousOrder, divisorData, config); - scope.Dispose(); - return order; - } + public static ulong Calculate( + ulong prime, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderSearchConfig config, + PrimeOrderHeuristicDevice device) + { + var scope = UsePow2Mode(device); + ulong order = CalculateInternal(prime, previousOrder, divisorData, config); + scope.Dispose(); + return order; + } - private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig 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 (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) - { - return gpuOrder; - } - - PartialFactorResult phiFactors = PartialFactor(phi, config); - - ulong result; - if (phiFactors.Factors is null) - { - result = FinishStrictlyCpu(prime, divisorData); - phiFactors.Dispose(); - return result; - } - - result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); - phiFactors.Dispose(); - return result; - } - - private static ulong RunHeuristicPipelineCpu( - ulong prime, - ulong? previousOrder, - in PrimeOrderSearchConfig 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, prime, phiFactors, divisorData); - - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config)) - { - return candidateOrder; - } - - if (config.Mode == PrimeOrderMode.Strict) - { - return FinishStrictlyCpu(prime, divisorData); - } - - if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) - { - return order; - } - - return candidateOrder; - } - - private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData) - { - return CalculateByFactorizationCpu(prime, divisorData); - } - - private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - for (int i = 0; i < length; i++) - { - ulong factor = factorSpan[i].Value; - ulong reduced = phi / factor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - - return true; - } - - 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, prime, divisorData)) - { - order = half; - } - } - - return order; - } - - private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - FactorEntry[]? tempArray = null; - try - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - tempArray = ArrayPool.Shared.Rent(length + 1); - Span buffer = tempArray; - factorSpan.CopyTo(buffer); - - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None); - if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime) - { - buffer[length] = new FactorEntry(factors.Cofactor, 1); - length++; - } - - buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value)); - - for (int i = 0; i < length; i++) - { - ulong primeFactor = buffer[i].Value; - int exponent = buffer[i].Exponent; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } - - ulong reduced = order / primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - order = reduced; - continue; - } - - break; - } - } - - return order; - } - finally - { - if (tempArray is not null) - { - ArrayPool.Shared.Return(tempArray, clearArray: false); - } - } - } - - private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) - { - if (order == 0UL) - { - return false; - } - - // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`. - // DebugLog("Verifying a^order ≡ 1 (mod p)"); - if (!Pow2EqualsOneCpu(order, prime, divisorData)) - { - return false; - } - - // DebugLog("Partial factoring order"); - - // TODO: Do we do partial factoring of order multiple times? - PartialFactorResult factorization = PartialFactor(order, config); - try - { - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL) - { - return false; - } - - // DebugLog("Cofactor > 1. Testing primality of cofactor"); - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); - if (!isPrime) - { - return false; - } - - // DebugLog("Adding cofactor as prime factor"); - PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); - factorization.Dispose(); - factorization = extended; - } - - ReadOnlySpan span = factorization.Factors; - // DebugLog("Verifying prime-power reductions"); - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - ulong primeFactor = span[i].Value; - ulong reduced = order; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - } - - return true; - } - finally - { - factorization.Dispose(); - } - } - - private static bool TryHeuristicFinishCpu( - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderSearchConfig config, - PartialFactorResult phiFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - return false; - } - - // TODO: Do we do partial factoring of order multiple times? - // DebugLog("Trying heuristic. Partial factoring order"); - PartialFactorResult orderFactors = PartialFactor(order, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL) - { - return false; - } - - var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor); - // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) - if (!isPrime) - { - return false; - } - - PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - orderFactors.Dispose(); - orderFactors = extended; - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; - List candidates = new(capacity); - FactorEntry[] factorArray = orderFactors.Factors!; - // DebugLog("Building candidates list"); - BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); - if (candidates.Count == 0) - { - return false; - } - - // DebugLog("Sorting candidates"); - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; - int powUsed = 0; - int candidateCount = candidates.Count; - bool allowGpuBatch = true; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - - // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); - int index = 0; - const int MaxGpuBatchSize = 256; - const int StackGpuBatchSize = 64; - Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; - while (index < candidateCount && powUsed < powBudget) - { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, 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 && IsGpuPow2Allowed) - { - if (batchSize <= StackGpuBatchSize) - { - Span localRemainders = stackGpuRemainders.Slice(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - gpuSuccess = true; - gpuStackRemainders = true; - } - } - else - { - gpuPool = ArrayPool.Shared.Rent(batchSize); - Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - pooledGpuRemainders = pooledRemainders; - gpuSuccess = true; - } - else - { - ArrayPool.Shared.Return(gpuPool, clearArray: false); - gpuPool = null; - } - } - - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) - { - allowGpuBatch = false; - } - } - - for (int i = 0; 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 = Pow2EqualsOneCpu(candidate, prime, divisorData); - } - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - if (gpuPool is not null) - { - ArrayPool.Shared.Return(gpuPool, clearArray: false); - } - - result = candidate; - return true; - } - - if (gpuPool is not null) - { - ArrayPool.Shared.Return(gpuPool, clearArray: false); - } - - index += batchSize; - } - - // DebugLog("No candidate confirmed"); - return false; - } - finally - { - orderFactors.Dispose(); - } - } - - private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) - { - ulong previous = previousOrder ?? 0UL; - int previousGroup = previousOrder.HasValue ? GetGroup(previousOrder.Value, prime) : 1; - bool hasPrevious = previousOrder.HasValue; - - candidates.Sort((x, y) => - { - CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious); - CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious); - int compare = keyX.Primary.CompareTo(keyY.Primary); - if (compare != 0) - { - return compare; - } - - compare = keyX.Secondary.CompareTo(keyY.Secondary); - if (compare != 0) - { - return compare; - } - - return keyX.Tertiary.CompareTo(keyY.Tertiary); - }); - } - - private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious) - { - int group = GetGroup(value, prime); - if (group == 0) - { - return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); - } - - bool isGe = !hasPrevious || value >= previous; - int primary = ComputePrimary(group, isGe, previousGroup); - long secondary; - long tertiary; - - if (group == 3) - { - secondary = -(long)value; - tertiary = -(long)value; - } - else - { - ulong reference = hasPrevious ? previous : 0UL; - ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; - secondary = (long)distance; - tertiary = (long)value; - } - - return new CandidateKey(primary, secondary, tertiary); - } - - 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 int GetGroup(ulong value, ulong prime) - { - ulong threshold1 = prime >> 3; - if (value <= threshold1) - { - return 1; - } - - ulong threshold2 = prime >> 2; - if (value <= threshold2) - { - return 2; - } - - ulong threshold3 = (prime * 3UL) >> 3; - if (value <= threshold3) - { - return 3; - } - - return 0; - } - - private static void BuildCandidates(ulong order, FactorEntry[] 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)); - BuildCandidatesRecursive(order, buffer, 0, 1UL, candidates, limit); - } - - private static void BuildCandidatesRecursive(ulong order, ReadOnlySpan factors, 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]; - ulong primeFactor = factor.Value; - ulong contribution = 1UL; - for (int exponent = 0; exponent <= factor.Exponent; exponent++) - { - ulong nextDivisor = divisorProduct * contribution; - if (nextDivisor > order) - { - break; - } - - BuildCandidatesRecursive(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 bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) - { - PartialFactorResult factorization = PartialFactor(candidate, config); - try - { - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL) - { - return false; - } - - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); - if (!isPrime) - { - return false; - } - - PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); - factorization.Dispose(); - factorization = extended; - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - ulong primeFactor = span[i].Value; - ulong reduced = candidate; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) - { - return false; - } - - powUsed++; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - } - - return true; - } - finally - { - factorization.Dispose(); - } - } - - private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) - { - if (IsGpuPow2Allowed) - { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); - if (status == GpuPow2ModStatus.Success) - { - return remainder == 1UL; - } - } - - return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; - } - - private static List AcquirePendingCompositeList(int capacityHint) - { - List? list = s_partialFactorPending; - if (list is null) - { - list = new List(capacityHint); - s_partialFactorPending = list; - return list; - } - - list.Clear(); - if (list.Capacity < capacityHint) - { - list.Capacity = capacityHint; - } - - return list; - } - - private static Stack AcquireCompositeStack(int capacityHint) - { - Stack? stack = s_partialFactorCompositeStack; - if (stack is null) - { - stack = new Stack(capacityHint); - s_partialFactorCompositeStack = stack; - return stack; - } - - stack.Clear(); - stack.EnsureCapacity(capacityHint); - return stack; - } - - private static Dictionary AcquireFactorCountDictionary(int capacityHint) - { - Dictionary? dictionary = s_partialFactorCounts; - if (dictionary is null) - { - dictionary = new Dictionary(capacityHint); - s_partialFactorCounts = dictionary; - return dictionary; - } - - dictionary.Clear(); - dictionary.EnsureCapacity(capacityHint); - return dictionary; - } - - private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - const int FactorSlotCount = GpuSmallPrimeFactorSlots; - Span primeSlots = stackalloc ulong[FactorSlotCount]; - Span exponentSlots = stackalloc int[FactorSlotCount]; - primeSlots.Clear(); - exponentSlots.Clear(); - - int factorCount = 0; - Dictionary? counts = null; - bool useDictionary = false; - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; - - bool gpuFactored = false; - if (IsGpuHeuristicDevice) - { - gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); - } - - List pending = AcquirePendingCompositeList(2); - Stack? compositeStack = null; - - try - { - if (!gpuFactored) - { - primeSlots.Clear(); - exponentSlots.Clear(); - if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) - { - useDictionary = true; - counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8)); - remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); - } - } - - if (remaining > 1UL) - { - pending.Add(remaining); - } - - if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) - { - Stopwatch stopwatch = Stopwatch.StartNew(); - long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks; - compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4)); - compositeStack.Push(remaining); - pending.Clear(); - - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - if (composite <= 1UL) - { - continue; - } - - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - if (isPrime) - { - AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - continue; - } - - if (stopwatch.ElapsedTicks > budgetTicks) - { - pending.Add(composite); - continue; - } - - if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor)) - { - pending.Add(composite); - continue; - } - - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - } - - ulong cofactor = 1UL; - int pendingCount = pending.Count; - for (int i = 0; i < pendingCount; i++) - { - ulong composite = pending[i]; - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - if (isPrime) - { - AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - } - else if (composite > 1UL) - { - cofactor = checked(cofactor * composite); - } - } - - if (useDictionary) - { - if ((counts is null || counts.Count == 0) && cofactor == value) - { - return PartialFactorResult.Rent(null, value, false, 0); - } - - if (counts is null || counts.Count == 0) - { - return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - } - - FactorEntry[] factors = ArrayPool.Shared.Rent(counts.Count); - int index = 0; - foreach (KeyValuePair entry in counts) - { - factors[index] = new FactorEntry(entry.Key, entry.Value); - index++; - } - - Span span = factors.AsSpan(0, index); - span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - - bool fullyFactored = cofactor == 1UL; - return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); - } - - int actualCount = 0; - for (int i = 0; i < factorCount; i++) - { - if (primeSlots[i] != 0UL && exponentSlots[i] != 0) - { - actualCount++; - } - } - - if (actualCount == 0) - { - if (cofactor == value) - { - return PartialFactorResult.Rent(null, value, false, 0); - } - - return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - } - - FactorEntry[] array = ArrayPool.Shared.Rent(actualCount); - int arrayIndex = 0; - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeSlots[i]; - int exponentValue = exponentSlots[i]; - if (primeValue == 0UL || exponentValue == 0) - { - continue; - } - - array[arrayIndex] = new FactorEntry(primeValue, exponentValue); - arrayIndex++; - } - - Span arraySpan = array.AsSpan(0, arrayIndex); - arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - bool fullyFactoredArray = cofactor == 1UL; - return PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); - } - finally - { - pending.Clear(); - counts?.Clear(); - compositeStack?.Clear(); - } - } - - private const int GpuSmallPrimeFactorSlots = 64; - - private static bool TryPopulateSmallPrimeFactorsCpu( - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining) - { - factorCount = 0; - remaining = value; - - if (value <= 1UL) - { - return true; - } - - int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); - if (capacity == 0) - { - return false; - } - - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; - ulong remainingLocal = value; - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - - for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > effectiveLimit) - { - break; - } - - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remainingLocal) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remainingLocal % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); - - if (factorCount >= capacity) - { - return false; - } - - primeTargets[factorCount] = primeValue; - exponentTargets[factorCount] = exponent; - factorCount++; - } - - remaining = remainingLocal; - return true; - } - - private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) - { - ulong 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; - } - - if (squares[i] > remaining) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == 0UL); - - counts[primeValue] = exponent; - } - - return remaining; - } - - private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) - { - Span primeBuffer = stackalloc ulong[GpuSmallPrimeFactorSlots]; - Span exponentBuffer = stackalloc int[GpuSmallPrimeFactorSlots]; - primeBuffer.Clear(); - exponentBuffer.Clear(); - - if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) - { - remaining = value; - return false; - } - - remaining = gpuRemaining; - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - if (primeValue == 0UL || exponent == 0) - { - continue; - } - - counts[primeValue] = exponent; - } - - return true; - } - private static bool TryPollardRho(ulong n, Stopwatch stopwatch, long budgetTicks, out ulong factor) - { - factor = 0UL; - if ((n & 1UL) == 0UL) - { - factor = 2UL; - return true; - } - - while (true) - { - if (stopwatch.ElapsedTicks > budgetTicks) - { - return false; - } - - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - if (stopwatch.ElapsedTicks > budgetTicks) - { - return false; - } - - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); - ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, n); - } - - if (d != n) - { - factor = d; - return true; - } - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) - { - UInt128 value = (UInt128)x * x + c; - return (ulong)(value % modulus); - } - - 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; - } - } - } - - 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; - } - } - - private static bool TryAppendFactor(Span primes, Span exponents, ref int count, ulong prime, int exponent) - { - if (exponent <= 0) - { - return true; - } - - for (int i = 0; i < count; i++) - { - if (primes[i] == prime) - { - exponents[i] += exponent; - return true; - } - } - - if (count >= primes.Length) - { - return false; - } - - primes[count] = prime; - exponents[count] = exponent; - count++; - return true; - } - - private static void AddFactorToCollector( - ref bool useDictionary, - ref Dictionary? counts, - Span primes, - Span exponents, - ref int count, - ulong prime, - int exponent) - { - if (exponent <= 0) - { - return; - } - - if (useDictionary) - { - AddFactor(counts!, prime, exponent); - return; - } - - if (TryAppendFactor(primes, exponents, ref count, prime, exponent)) - { - return; - } - - counts ??= AcquireFactorCountDictionary(Math.Max(count, 8)); - CopyFactorsToDictionary(primes, exponents, count, counts); - count = 0; - useDictionary = true; - AddFactor(counts, prime, exponent); - } - - private static void CopyFactorsToDictionary( - ReadOnlySpan primes, - ReadOnlySpan exponents, - int count, - Dictionary target) - { - for (int i = 0; i < count; i++) - { - ulong prime = primes[i]; - int exponent = exponents[i]; - if (prime == 0UL || exponent == 0) - { - continue; - } - - target[prime] = exponent; - } - } - - private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) - { - ulong phi = prime - 1UL; - Dictionary counts = new(capacity: 8); - FactorCompletely(phi, counts); - if (counts.Count == 0) - { - return phi; - } - - List> entries = new(counts); - entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); - - ulong order = phi; - int entryCount = entries.Count; - for (int i = 0; i < entryCount; i++) - { - ulong primeFactor = entries[i].Key; - int exponent = entries[i].Value; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } - - ulong candidate = order / primeFactor; - if (Pow2EqualsOneCpu(candidate, prime, divisorData)) - { - order = candidate; - continue; - } - - break; - } - } - - return order; - } - - private static void FactorCompletely(ulong value, Dictionary counts) - { - if (value <= 1UL) - { - return; - } - - if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value)) - { - AddFactor(counts, value, 1); - return; - } - - ulong factor = PollardRhoStrict(value); - FactorCompletely(factor, counts); - FactorCompletely(value / factor, counts); - } - - private static ulong PollardRhoStrict(ulong n) - { - if ((n & 1UL) == 0UL) - { - return 2UL; - } - - while (true) - { - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); - ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, n); - } - - if (d != n) - { - return d; - } - } - } - - private readonly struct CandidateKey(int primary, long secondary, long tertiary) - { - public readonly int Primary = primary; - public readonly long Secondary = secondary; - public readonly long Tertiary = tertiary; - } + if (prime <= 3UL) + { + return prime == 3UL ? 2UL : 1UL; + } + + ulong phi = prime - 1UL; + + if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) + { + return gpuOrder; + } + + PartialFactorResult phiFactors = PartialFactor(phi, config); + + ulong result; + if (phiFactors.Factors is null) + { + result = FinishStrictlyCpu(prime, divisorData); + phiFactors.Dispose(); + return result; + } + + result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); + phiFactors.Dispose(); + return result; + } + + private static ulong RunHeuristicPipelineCpu( + ulong prime, + ulong? previousOrder, + in PrimeOrderSearchConfig 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, prime, phiFactors, divisorData); + + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config)) + { + return candidateOrder; + } + + if (config.Mode == PrimeOrderMode.Strict) + { + return FinishStrictlyCpu(prime, divisorData); + } + + if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) + { + return order; + } + + return candidateOrder; + } + + private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData) + { + return CalculateByFactorizationCpu(prime, divisorData); + } + + private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + for (int i = 0; i < length; i++) + { + ulong factor = factorSpan[i].Value; + ulong reduced = phi / factor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; + } + } + + return true; + } + + 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, prime, divisorData)) + { + order = half; + } + } + + return order; + } + + private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + FactorEntry[] tempArray = pool.Rent(length + 1); + + Span buffer = tempArray; + factorSpan.CopyTo(buffer); + + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None); + if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime) + { + buffer[length] = new FactorEntry(factors.Cofactor, 1); + length++; + } + + buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value)); + + for (int i = 0; i < length; i++) + { + ulong primeFactor = buffer[i].Value; + int exponent = buffer[i].Exponent; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong reduced = order / primeFactor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + order = reduced; + continue; + } + + break; + } + } + + pool.Return(tempArray, clearArray: false); + return order; + } + + private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + { + if (order == 0UL) + { + return false; + } + + // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`. + // DebugLog("Verifying a^order ≡ 1 (mod p)"); + if (!Pow2EqualsOneCpu(order, prime, divisorData)) + { + return false; + } + + // DebugLog("Partial factoring order"); + + // TODO: Do we do partial factoring of order multiple times? + PartialFactorResult factorization = PartialFactor(order, config); + try + { + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) + { + return false; + } + + // DebugLog("Cofactor > 1. Testing primality of cofactor"); + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); + if (!isPrime) + { + return false; + } + + // DebugLog("Adding cofactor as prime factor"); + PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.Dispose(); + factorization = extended; + } + + ReadOnlySpan span = factorization.Factors; + // DebugLog("Verifying prime-power reductions"); + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + ulong primeFactor = span[i].Value; + ulong reduced = order; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; + } + } + } + + return true; + } + finally + { + factorization.Dispose(); + } + } + + private static bool TryHeuristicFinishCpu( + ulong prime, + ulong order, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderSearchConfig config, + PartialFactorResult phiFactors, + out ulong result) + { + result = 0UL; + if (order <= 1UL) + { + return false; + } + + // TODO: Do we do partial factoring of order multiple times? + // DebugLog("Trying heuristic. Partial factoring order"); + PartialFactorResult orderFactors = PartialFactor(order, config); + try + { + if (orderFactors.Factors is null) + { + return false; + } + + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= 1UL) + { + return false; + } + + var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor); + // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) + if (!isPrime) + { + return false; + } + + PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + orderFactors.Dispose(); + orderFactors = extended; + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; + List candidates = new(capacity); + FactorEntry[] factorArray = orderFactors.Factors!; + // DebugLog("Building candidates list"); + BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) + { + return false; + } + + // DebugLog("Sorting candidates"); + SortCandidates(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; + int powUsed = 0; + int candidateCount = candidates.Count; + bool allowGpuBatch = true; + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + + // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); + int index = 0; + const int MaxGpuBatchSize = 256; + const int StackGpuBatchSize = 64; + Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; + 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, 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 && IsGpuPow2Allowed) + { + if (batchSize <= StackGpuBatchSize) + { + Span localRemainders = stackGpuRemainders.Slice(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(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(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; + } + } + + for (int i = 0; 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 = Pow2EqualsOneCpu(candidate, prime, divisorData); + } + if (!equalsOne) + { + continue; + } + + if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + continue; + } + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + result = candidate; + return true; + } + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + index += batchSize; + } + + // DebugLog("No candidate confirmed"); + return false; + } + finally + { + orderFactors.Dispose(); + } + } + + private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) + { + ulong previous = previousOrder ?? 0UL; + int previousGroup = previousOrder.HasValue ? GetGroup(previousOrder.Value, prime) : 1; + bool hasPrevious = previousOrder.HasValue; + + candidates.Sort((x, y) => + { + CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious); + CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious); + int compare = keyX.Primary.CompareTo(keyY.Primary); + if (compare != 0) + { + return compare; + } + + compare = keyX.Secondary.CompareTo(keyY.Secondary); + if (compare != 0) + { + return compare; + } + + return keyX.Tertiary.CompareTo(keyY.Tertiary); + }); + } + + private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious) + { + int group = GetGroup(value, prime); + if (group == 0) + { + return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); + } + + bool isGe = !hasPrevious || value >= previous; + int primary = ComputePrimary(group, isGe, previousGroup); + long secondary; + long tertiary; + + if (group == 3) + { + secondary = -(long)value; + tertiary = -(long)value; + } + else + { + ulong reference = hasPrevious ? previous : 0UL; + ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; + secondary = (long)distance; + tertiary = (long)value; + } + + return new CandidateKey(primary, secondary, tertiary); + } + + 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 int GetGroup(ulong value, ulong prime) + { + ulong threshold1 = prime >> 3; + if (value <= threshold1) + { + return 1; + } + + ulong threshold2 = prime >> 2; + if (value <= threshold2) + { + return 2; + } + + ulong threshold3 = (prime * 3UL) >> 3; + if (value <= threshold3) + { + return 3; + } + + return 0; + } + + private static void BuildCandidates(ulong order, FactorEntry[] 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)); + BuildCandidatesRecursive(order, buffer, 0, 1UL, candidates, limit); + } + + private static void BuildCandidatesRecursive(ulong order, in ReadOnlySpan factors, 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]; + ulong primeFactor = factor.Value; + ulong contribution = 1UL; + for (int exponent = 0; exponent <= factor.Exponent; exponent++) + { + ulong nextDivisor = divisorProduct * contribution; + if (nextDivisor > order) + { + break; + } + + BuildCandidatesRecursive(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 bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) + { + PartialFactorResult factorization = PartialFactor(candidate, config); + try + { + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) + { + return false; + } + + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); + if (!isPrime) + { + return false; + } + + PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.Dispose(); + factorization = extended; + } + + ReadOnlySpan span = factorization.Factors; + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + ulong primeFactor = span[i].Value; + ulong reduced = candidate; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (powUsed >= powBudget && powBudget > 0) + { + return false; + } + + powUsed++; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; + } + } + } + + return true; + } + finally + { + factorization.Dispose(); + } + } + + private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) + { + if (IsGpuPow2Allowed) + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); + if (status == GpuPow2ModStatus.Success) + { + return remainder == 1UL; + } + } + + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; + } + + private static List AcquirePendingCompositeList(int capacityHint) + { + List? list = s_partialFactorPending; + if (list is null) + { + list = new List(capacityHint); + s_partialFactorPending = list; + return list; + } + + list.Clear(); + if (list.Capacity < capacityHint) + { + list.Capacity = capacityHint; + } + + return list; + } + + private static Stack AcquireCompositeStack(int capacityHint) + { + Stack? stack = s_partialFactorCompositeStack; + if (stack is null) + { + stack = new Stack(capacityHint); + s_partialFactorCompositeStack = stack; + return stack; + } + + stack.Clear(); + stack.EnsureCapacity(capacityHint); + return stack; + } + + private static Dictionary AcquireFactorCountDictionary(int capacityHint) + { + Dictionary? dictionary = s_partialFactorCounts; + if (dictionary is null) + { + dictionary = new Dictionary(capacityHint); + s_partialFactorCounts = dictionary; + return dictionary; + } + + dictionary.Clear(); + dictionary.EnsureCapacity(capacityHint); + return dictionary; + } + + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + const int FactorSlotCount = GpuSmallPrimeFactorSlots; + Span primeSlots = stackalloc ulong[FactorSlotCount]; + Span exponentSlots = stackalloc int[FactorSlotCount]; + primeSlots.Clear(); + exponentSlots.Clear(); + + int factorCount = 0; + Dictionary? counts = null; + bool useDictionary = false; + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong remaining = value; + + bool gpuFactored = false; + if (IsGpuHeuristicDevice) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); + } + + List pending = AcquirePendingCompositeList(2); + Stack? compositeStack = null; + + try + { + if (!gpuFactored) + { + primeSlots.Clear(); + exponentSlots.Clear(); + if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) + { + useDictionary = true; + counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8)); + remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + } + } + + if (remaining > 1UL) + { + pending.Add(remaining); + } + + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + Stopwatch stopwatch = Stopwatch.StartNew(); + long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks; + compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4)); + compositeStack.Push(remaining); + pending.Clear(); + + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + if (composite <= 1UL) + { + continue; + } + + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + continue; + } + + if (stopwatch.ElapsedTicks > budgetTicks) + { + pending.Add(composite); + continue; + } + + if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor)) + { + pending.Add(composite); + continue; + } + + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + } + + ulong cofactor = 1UL; + int pendingCount = pending.Count; + for (int i = 0; i < pendingCount; i++) + { + ulong composite = pending[i]; + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + } + else if (composite > 1UL) + { + cofactor = checked(cofactor * composite); + } + } + + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + if (useDictionary) + { + if ((counts is null || counts.Count == 0) && cofactor == value) + { + return PartialFactorResult.Rent(null, value, false, 0); + } + + if (counts is null || counts.Count == 0) + { + return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + } + + FactorEntry[] factors = pool.Rent(counts.Count); + int index = 0; + foreach (KeyValuePair entry in counts) + { + factors[index] = new FactorEntry(entry.Key, entry.Value); + index++; + } + + Span span = factors.AsSpan(0, index); + span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + + bool fullyFactored = cofactor == 1UL; + return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); + } + + int actualCount = 0; + for (int i = 0; i < factorCount; i++) + { + if (primeSlots[i] != 0UL && exponentSlots[i] != 0) + { + actualCount++; + } + } + + if (actualCount == 0) + { + if (cofactor == value) + { + return PartialFactorResult.Rent(null, value, false, 0); + } + + return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + } + + FactorEntry[] array = pool.Rent(actualCount); + int arrayIndex = 0; + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeSlots[i]; + int exponentValue = exponentSlots[i]; + if (primeValue == 0UL || exponentValue == 0) + { + continue; + } + + array[arrayIndex] = new FactorEntry(primeValue, exponentValue); + arrayIndex++; + } + + Span arraySpan = array.AsSpan(0, arrayIndex); + arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + bool fullyFactoredArray = cofactor == 1UL; + return PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); + } + finally + { + pending.Clear(); + counts?.Clear(); + compositeStack?.Clear(); + } + } + + private const int GpuSmallPrimeFactorSlots = 64; + + private static bool TryPopulateSmallPrimeFactorsCpu( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) + { + factorCount = 0; + remaining = value; + + if (value <= 1UL) + { + return true; + } + + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + if (capacity == 0) + { + return false; + } + + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + ulong remainingLocal = value; + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remainingLocal) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remainingLocal % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); + + if (factorCount >= capacity) + { + return false; + } + + primeTargets[factorCount] = primeValue; + exponentTargets[factorCount] = exponent; + factorCount++; + } + + remaining = remainingLocal; + return true; + } + + private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) + { + ulong 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; + } + + if (squares[i] > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + counts[primeValue] = exponent; + } + + return remaining; + } + + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) + { + Span primeBuffer = stackalloc ulong[GpuSmallPrimeFactorSlots]; + Span exponentBuffer = stackalloc int[GpuSmallPrimeFactorSlots]; + primeBuffer.Clear(); + exponentBuffer.Clear(); + + if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) + { + remaining = value; + return false; + } + + remaining = gpuRemaining; + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + if (primeValue == 0UL || exponent == 0) + { + continue; + } + + counts[primeValue] = exponent; + } + + return true; + } + private static bool TryPollardRho(ulong n, Stopwatch stopwatch, long budgetTicks, out ulong factor) + { + factor = 0UL; + if ((n & 1UL) == 0UL) + { + factor = 2UL; + return true; + } + + while (true) + { + if (stopwatch.ElapsedTicks > budgetTicks) + { + return false; + } + + ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + if (stopwatch.ElapsedTicks > budgetTicks) + { + return false; + } + + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomial(y, c, n); + y = AdvancePolynomial(y, c, n); + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, n); + } + + if (d != n) + { + factor = d; + return true; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) + { + UInt128 value = (UInt128)x * x + c; + return (ulong)(value % modulus); + } + + 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; + } + } + } + + 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; + } + } + + private static bool TryAppendFactor(Span primes, Span exponents, ref int count, ulong prime, int exponent) + { + if (exponent <= 0) + { + return true; + } + + for (int i = 0; i < count; i++) + { + if (primes[i] == prime) + { + exponents[i] += exponent; + return true; + } + } + + if (count >= primes.Length) + { + return false; + } + + primes[count] = prime; + exponents[count] = exponent; + count++; + return true; + } + + private static void AddFactorToCollector( + ref bool useDictionary, + ref Dictionary? counts, + Span primes, + Span exponents, + ref int count, + ulong prime, + int exponent) + { + if (exponent <= 0) + { + return; + } + + if (useDictionary) + { + AddFactor(counts!, prime, exponent); + return; + } + + if (TryAppendFactor(primes, exponents, ref count, prime, exponent)) + { + return; + } + + counts ??= AcquireFactorCountDictionary(Math.Max(count, 8)); + CopyFactorsToDictionary(primes, exponents, count, counts); + count = 0; + useDictionary = true; + AddFactor(counts, prime, exponent); + } + + private static void CopyFactorsToDictionary( + in ReadOnlySpan primes, + in ReadOnlySpan exponents, + int count, + Dictionary target) + { + for (int i = 0; i < count; i++) + { + ulong prime = primes[i]; + int exponent = exponents[i]; + if (prime == 0UL || exponent == 0) + { + continue; + } + + target[prime] = exponent; + } + } + + private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) + { + ulong phi = prime - 1UL; + Dictionary counts = new(capacity: 8); + FactorCompletely(phi, counts); + if (counts.Count == 0) + { + return phi; + } + + List> entries = [.. counts]; + entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); + + ulong order = phi; + int entryCount = entries.Count; + for (int i = 0; i < entryCount; i++) + { + ulong primeFactor = entries[i].Key; + int exponent = entries[i].Value; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong candidate = order / primeFactor; + if (Pow2EqualsOneCpu(candidate, prime, divisorData)) + { + order = candidate; + continue; + } + + break; + } + } + + return order; + } + + private static void FactorCompletely(ulong value, Dictionary counts) + { + if (value <= 1UL) + { + return; + } + + if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value)) + { + AddFactor(counts, value, 1); + return; + } + + ulong factor = PollardRhoStrict(value); + FactorCompletely(factor, counts); + FactorCompletely(value / factor, counts); + } + + private static ulong PollardRhoStrict(ulong n) + { + if ((n & 1UL) == 0UL) + { + return 2UL; + } + + while (true) + { + ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomial(y, c, n); + y = AdvancePolynomial(y, c, n); + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, n); + } + + if (d != n) + { + return d; + } + } + } + + private readonly struct CandidateKey(int primary, long secondary, long tertiary) + { + public readonly int Primary = primary; + public readonly long Secondary = secondary; + public readonly long Tertiary = tertiary; + } } diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs new file mode 100644 index 0000000..cd02a0d --- /dev/null +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -0,0 +1,41 @@ +using System.Buffers; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core +{ + public static class ThreadStaticPools + { + [ThreadStatic] + private static ArrayPool? _factorEntryPool; + public static ArrayPool FactorEntryPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _factorEntryPool ??= ArrayPool.Create(); + } + } + + [ThreadStatic] + private static ArrayPool? _factorEntry128Pool; + public static ArrayPool FactorEntry128Pool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _factorEntry128Pool ??= ArrayPool.Create(); + } + } + + [ThreadStatic] + private static ArrayPool? _ulongPool; + public static ArrayPool UlongPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _ulongPool ??= ArrayPool.Create(); + } + } + } +} \ No newline at end of file