From 51730bcba323b47b8d642f380fb25a9124c909ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 15 Oct 2025 23:43:27 +0200 Subject: [PATCH 01/14] Reduce GC pressure in divisor checking pipeline --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 220 +++++++++++--------- PerfectNumbers.Core/ULongExtensions.cs | 105 ++++++---- 2 files changed, 183 insertions(+), 142 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 39d710c..6cc4f33 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -1,4 +1,5 @@ using System.Buffers; +using System.Collections.Generic; using System.Buffers.Binary; using System.Diagnostics; using System.Numerics; @@ -56,6 +57,9 @@ internal static partial class PrimeOrderCalculator [ThreadStatic] private static Dictionary? s_partialFactorCounts; + [ThreadStatic] + private static List? s_heuristicCandidateList; + private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; private readonly struct Pow2ModeScope @@ -402,124 +406,131 @@ internal static partial class PrimeOrderCalculator } 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) + List candidates = AcquireHeuristicCandidateList(capacity); + try { - 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) + FactorEntry[] factorArray = orderFactors.Factors!; + // DebugLog("Building candidates list"); + BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) { - break; + return false; } - - 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) + + // 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) { - if (batchSize <= StackGpuBatchSize) + int remaining = candidateCount - index; + int budgetRemaining = powBudget - powUsed; + int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize)); + if (batchSize <= 0) { - Span localRemainders = stackGpuRemainders.Slice(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - gpuSuccess = true; - gpuStackRemainders = true; - } + break; } - else + + 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) { - gpuPool = pool.Rent(batchSize); - Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) + if (batchSize <= StackGpuBatchSize) { - pooledGpuRemainders = pooledRemainders; - gpuSuccess = true; + Span localRemainders = stackGpuRemainders.Slice(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + gpuSuccess = true; + gpuStackRemainders = true; + } } else { - pool.Return(gpuPool, clearArray: false); - gpuPool = null; + 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; } } - - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) + + for (int i = 0; i < batchSize && powUsed < powBudget; i++) { - allowGpuBatch = false; + 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; } - } - - 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; + + index += batchSize; } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - index += batchSize; + + } + finally + { + candidates.Clear(); } - // DebugLog("No candidate confirmed"); return false; } @@ -823,6 +834,21 @@ internal static partial class PrimeOrderCalculator return dictionary; } + private static List AcquireHeuristicCandidateList(int capacityHint) + { + List? list = s_heuristicCandidateList; + if (list is null) + { + list = new List(capacityHint); + s_heuristicCandidateList = list; + return list; + } + + list.Clear(); + list.EnsureCapacity(capacityHint); + return list; + } + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) { if (value <= 1UL) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index fa0d030..9edc591 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -346,48 +346,67 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount); + + const int StackOddPowerThreshold = PerfectNumberConstants.PooledArrayThreshold; + Span oddPowersSpan = oddPowerCount <= StackOddPowerThreshold + ? stackalloc ulong[StackOddPowerThreshold] + : default; + + Span oddPowers; + ulong[]? oddPowersArray = null; + + if (!oddPowersSpan.IsEmpty) + { + oddPowers = oddPowersSpan.Slice(0, oddPowerCount); + } + else + { + oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); + } + + InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers); int index = bitLength - 1; - while (index >= 0) - { - if (((exponent >> index) & 1UL) == 0UL) - { - result = result.MontgomeryMultiply(result, modulus, nPrime); - index--; - continue; - } + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + index--; + continue; + } - int windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } - while (((exponent >> windowStart) & 1UL) == 0UL) - { - windowStart++; - } + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } - int windowLength = index - windowStart + 1; - for (int square = 0; square < windowLength; square++) - { - result = result.MontgomeryMultiply(result, modulus, nPrime); - } + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + } - ulong mask = (1UL << windowLength) - 1UL; - ulong windowValue = (exponent >> windowStart) & mask; - int tableIndex = (int)((windowValue - 1UL) >> 1); - ulong multiplier = oddPowers[tableIndex]; - result = result.MontgomeryMultiply(multiplier, modulus, nPrime); + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MontgomeryMultiply(multiplier, modulus, nPrime); - index = windowStart - 1; - } - - if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold) - { - ArrayPool.Shared.Return(oddPowers); - } + index = windowStart - 1; + } + + if (oddPowersArray is not null) + { + ThreadStaticPools.UlongPool.Return(oddPowersArray, clearArray: false); + } if (keepMontgomery) { @@ -533,26 +552,22 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span oddPowers) { - ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold ? new ulong[oddPowerCount] : ArrayPool.Shared.Rent(oddPowerCount); - oddPowers[0] = divisor.MontgomeryTwo; - if (oddPowerCount == 1) + if (oddPowers.Length == 1) { - return oddPowers; + return; } - ulong previous, - square = divisor.MontgomeryTwoSquared; + ulong previous; + ulong square = divisor.MontgomeryTwoSquared; - for (int i = 1; i < oddPowerCount; i++) + for (int i = 1; i < oddPowers.Length; i++) { previous = oddPowers[i - 1]; oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); } - - return oddPowers; } [MethodImpl(MethodImplOptions.AggressiveInlining)] From 0a2350244139fdeee01b6538727eae081157b179 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 15 Oct 2025 23:57:46 +0200 Subject: [PATCH 02/14] Clear heuristic candidate list without try-finally --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 202 ++++++++++---------- 1 file changed, 99 insertions(+), 103 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 6cc4f33..ec6f006 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -407,131 +407,127 @@ internal static partial class PrimeOrderCalculator int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; List candidates = AcquireHeuristicCandidateList(capacity); - try + FactorEntry[] factorArray = orderFactors.Factors!; + // DebugLog("Building candidates list"); + BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) { - FactorEntry[] factorArray = orderFactors.Factors!; - // DebugLog("Building candidates list"); - BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); - if (candidates.Count == 0) + candidates.Clear(); + 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) { - return false; + break; } - - // 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) + + 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) { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize)); - if (batchSize <= 0) + if (batchSize <= StackGpuBatchSize) { - 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) { - 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; + gpuSuccess = true; + gpuStackRemainders = true; } } - - for (int i = 0; i < batchSize && powUsed < powBudget; i++) + else { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne; - if (gpuSuccess) + gpuPool = pool.Rent(batchSize); + Span pooledRemainders = gpuPool.AsSpan(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) { - ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; - equalsOne = remainderValue == 1UL; + pooledGpuRemainders = pooledRemainders; + gpuSuccess = true; } 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); + gpuPool = null; } - - result = candidate; - return true; } - + + 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); } - - index += batchSize; + + candidates.Clear(); + result = candidate; + return true; } - + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + index += batchSize; } - finally - { - candidates.Clear(); - } - // DebugLog("No candidate confirmed"); + + candidates.Clear(); + // DebugLog("No candidate confirmed"); return false; } finally From 3e7bffc8127d6006f72eb9e71223a8e75661d6fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 00:37:01 +0200 Subject: [PATCH 03/14] Refactor heuristic pooling to cut GC pressure --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 523 +++++++++----------- PerfectNumbers.Core/ThreadStaticPools.cs | 124 +++++ PerfectNumbers.Core/ULongExtensions.cs | 24 +- 3 files changed, 356 insertions(+), 315 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index ec6f006..8fddc0b 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -48,17 +48,9 @@ internal static partial class PrimeOrderCalculator [ThreadStatic] private static bool s_debugLoggingEnabled; - [ThreadStatic] - private static List? s_partialFactorPending; - [ThreadStatic] - private static Stack? s_partialFactorCompositeStack; - [ThreadStatic] - private static Dictionary? s_partialFactorCounts; - [ThreadStatic] - private static List? s_heuristicCandidateList; private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; @@ -406,13 +398,13 @@ internal static partial class PrimeOrderCalculator } int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; - List candidates = AcquireHeuristicCandidateList(capacity); + List candidates = ThreadStaticPools.RentUlongList(capacity); FactorEntry[] factorArray = orderFactors.Factors!; // DebugLog("Building candidates list"); BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); if (candidates.Count == 0) { - candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); return false; } @@ -513,7 +505,7 @@ internal static partial class PrimeOrderCalculator pool.Return(gpuPool, clearArray: false); } - candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); result = candidate; return true; } @@ -526,7 +518,7 @@ internal static partial class PrimeOrderCalculator index += batchSize; } - candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); // DebugLog("No candidate confirmed"); return false; } @@ -781,324 +773,271 @@ internal static partial class PrimeOrderCalculator 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 List AcquireHeuristicCandidateList(int capacityHint) - { - List? list = s_heuristicCandidateList; - if (list is null) - { - list = new List(capacityHint); - s_heuristicCandidateList = list; - return list; - } - - list.Clear(); - list.EnsureCapacity(capacityHint); - return list; - } - private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } - const int FactorSlotCount = GpuSmallPrimeFactorSlots; - Span primeSlots = stackalloc ulong[FactorSlotCount]; - Span exponentSlots = stackalloc int[FactorSlotCount]; - primeSlots.Clear(); - exponentSlots.Clear(); + 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; + int factorCount = 0; + Dictionary? counts = null; + bool useDictionary = false; - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; + 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 _); - } + bool gpuFactored = false; + if (IsGpuHeuristicDevice) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); + } - List pending = AcquirePendingCompositeList(2); - Stack? compositeStack = null; + List pending = ThreadStaticPools.RentUlongList(2); + Stack? compositeStack = null; + PartialFactorResult result; - 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 (!gpuFactored) + { + primeSlots.Clear(); + exponentSlots.Clear(); + if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) + { + useDictionary = true; + counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(factorCount, 8)); + remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + } + } - if (remaining > 1UL) - { - pending.Add(remaining); - } + if (remaining > 1UL) + { + pending.Add(remaining); + } - if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4)); - compositeStack.Push(remaining); - pending.Clear(); + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count * 2, 4)); + compositeStack.Push(remaining); + pending.Clear(); - long timestamp = 0L; // reused for deadline checks. - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - if (composite <= 1UL) - { - continue; - } + long timestamp = 0L; + 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; - } + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + continue; + } - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - pending.Add(composite); - continue; - } + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + pending.Add(composite); + continue; + } - if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(composite); - continue; - } + if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(composite); + continue; + } - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - } + 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); - } - } + 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); - } + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + if (useDictionary) + { + if ((counts is null || counts.Count == 0) && cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; + } - if (counts is null || counts.Count == 0) - { - return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - } + if (counts is null || counts.Count == 0) + { + result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + goto ReturnResult; + } - FactorEntry[] factors = pool.Rent(counts.Count); - int index = 0; - foreach (KeyValuePair entry in counts) - { - factors[index] = new FactorEntry(entry.Key, entry.Value); - index++; - } + 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)); + 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); - } + bool fullyFactored = cofactor == 1UL; + result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); + goto ReturnResult; + } - int actualCount = 0; - for (int i = 0; i < factorCount; i++) - { - if (primeSlots[i] != 0UL && exponentSlots[i] != 0) - { - actualCount++; - } - } + 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); - } + if (actualCount == 0) + { + if (cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; + } - return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - } + result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + goto ReturnResult; + } - 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; - } + 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++; - } + 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(); - } - } + Span arraySpan = array.AsSpan(0, arrayIndex); + arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + bool fullyFactoredArray = cofactor == 1UL; + result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); - private const int GpuSmallPrimeFactorSlots = 64; + ReturnResult: + if (counts is not null) + { + ThreadStaticPools.ReturnUlongIntDictionary(counts); + } - private static bool TryPopulateSmallPrimeFactorsCpu( - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining) - { - factorCount = 0; - remaining = value; + ThreadStaticPools.ReturnUlongList(pending); + if (compositeStack is not null) + { + ThreadStaticPools.ReturnUlongStack(compositeStack); + } - if (value <= 1UL) - { - return true; - } + return result; + } - int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); - if (capacity == 0) - { - return false; - } + private const int GpuSmallPrimeFactorSlots = 64; - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; - ulong remainingLocal = value; - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + private static bool TryPopulateSmallPrimeFactorsCpu( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) + { + factorCount = 0; + remaining = value; - for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > effectiveLimit) - { - break; - } + if (value <= 1UL) + { + return true; + } - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remainingLocal) - { - break; - } + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + if (capacity == 0) + { + return false; + } - ulong primeValue = primeCandidate; - if ((remainingLocal % primeValue) != 0UL) - { - continue; - } + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + ulong remainingLocal = value; + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } - if (factorCount >= capacity) - { - return false; - } + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remainingLocal) + { + break; + } - primeTargets[factorCount] = primeValue; - exponentTargets[factorCount] = exponent; - factorCount++; - } + ulong primeValue = primeCandidate; + if ((remainingLocal % primeValue) != 0UL) + { + continue; + } - remaining = remainingLocal; - return true; - } + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); - private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) + 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; @@ -1340,7 +1279,7 @@ internal static partial class PrimeOrderCalculator return; } - counts ??= AcquireFactorCountDictionary(Math.Max(count, 8)); + counts ??= ThreadStaticPools.RentUlongIntDictionary(Math.Max(count, 8)); CopyFactorsToDictionary(primes, exponents, count, counts); count = 0; useDictionary = true; diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 313e30b..9ba6060 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -79,6 +79,130 @@ namespace PerfectNumbers.Core } } + + [ThreadStatic] + private static List>? _ulongListPool; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static List RentUlongList(int capacityHint) + { + List>? pool = _ulongListPool; + if (pool is not null && pool.Count > 0) + { + int lastIndex = pool.Count - 1; + List list = pool[lastIndex]; + pool.RemoveAt(lastIndex); + if (list.Count > 0) + { + list.Clear(); + } + + if (list.Capacity < capacityHint) + { + list.Capacity = capacityHint; + } + + return list; + } + + return new List(capacityHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void ReturnUlongList(List list) + { + List>? pool = _ulongListPool; + if (pool is null) + { + pool = new List>(4); + _ulongListPool = pool; + } + + pool.Add(list); + } + + [ThreadStatic] + private static List>? _ulongStackPool; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Stack RentUlongStack(int capacityHint) + { + List>? pool = _ulongStackPool; + if (pool is not null && pool.Count > 0) + { + int lastIndex = pool.Count - 1; + Stack stack = pool[lastIndex]; + pool.RemoveAt(lastIndex); + if (stack.Count > 0) + { + stack.Clear(); + } + + if (stack.Count < capacityHint) + { + stack.EnsureCapacity(capacityHint); + } + + return stack; + } + + return new Stack(capacityHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void ReturnUlongStack(Stack stack) + { + List>? pool = _ulongStackPool; + if (pool is null) + { + pool = new List>(4); + _ulongStackPool = pool; + } + + pool.Add(stack); + } + + [ThreadStatic] + private static List>? _ulongIntDictionaryPool; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Dictionary RentUlongIntDictionary(int capacityHint) + { + List>? pool = _ulongIntDictionaryPool; + if (pool is not null && pool.Count > 0) + { + int lastIndex = pool.Count - 1; + Dictionary dictionary = pool[lastIndex]; + pool.RemoveAt(lastIndex); + if (dictionary.Count > 0) + { + dictionary.Clear(); + } + + if (dictionary.Count < capacityHint) + { + dictionary.EnsureCapacity(capacityHint); + } + + return dictionary; + } + + return new Dictionary(capacityHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void ReturnUlongIntDictionary(Dictionary dictionary) + { + List>? pool = _ulongIntDictionaryPool; + if (pool is null) + { + pool = new List>(4); + _ulongIntDictionaryPool = pool; + } + + pool.Add(dictionary); + } + [ThreadStatic] private static Dictionary? _mersenneFactorCacheDictionary; diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 9edc591..0271021 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -1,5 +1,4 @@ using System; -using System.Buffers; using System.Numerics; using System.Runtime.CompilerServices; using ILGPU.Algorithms; @@ -347,24 +346,7 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - const int StackOddPowerThreshold = PerfectNumberConstants.PooledArrayThreshold; - Span oddPowersSpan = oddPowerCount <= StackOddPowerThreshold - ? stackalloc ulong[StackOddPowerThreshold] - : default; - - Span oddPowers; - ulong[]? oddPowersArray = null; - - if (!oddPowersSpan.IsEmpty) - { - oddPowers = oddPowersSpan.Slice(0, oddPowerCount); - } - else - { - oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); - oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); - } - + Span oddPowers = stackalloc ulong[oddPowerCount]; InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers); int index = bitLength - 1; @@ -403,10 +385,6 @@ public static class ULongExtensions index = windowStart - 1; } - if (oddPowersArray is not null) - { - ThreadStaticPools.UlongPool.Return(oddPowersArray, clearArray: false); - } if (keepMontgomery) { From 9c474fc241ac5be8cf5e09f198af913d49640f76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 00:59:00 +0200 Subject: [PATCH 04/14] Use stack buffers for small composite collections --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 191 +++++++++++++++++++- 1 file changed, 182 insertions(+), 9 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 8fddc0b..da15d5c 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -799,8 +799,8 @@ internal static partial class PrimeOrderCalculator gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); } - List pending = ThreadStaticPools.RentUlongList(2); - Stack? compositeStack = null; + Span pendingStorage = stackalloc ulong[4]; + PendingBuffer pending = new PendingBuffer(pendingStorage); PartialFactorResult result; if (!gpuFactored) @@ -822,8 +822,10 @@ internal static partial class PrimeOrderCalculator if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) { + Span compositeStorage = stackalloc ulong[8]; + CompositeStack compositeStack = new CompositeStack(compositeStorage); + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count * 2, 4)); compositeStack.Push(remaining); pending.Clear(); @@ -860,13 +862,16 @@ internal static partial class PrimeOrderCalculator compositeStack.Push(factor); compositeStack.Push(quotient); } + + compositeStack.Dispose(); } ulong cofactor = 1UL; - int pendingCount = pending.Count; + ReadOnlySpan pendingSpan = pending.AsSpan(); + int pendingCount = pendingSpan.Length; for (int i = 0; i < pendingCount; i++) { - ulong composite = pending[i]; + ulong composite = pendingSpan[i]; bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) { @@ -956,13 +961,181 @@ internal static partial class PrimeOrderCalculator ThreadStaticPools.ReturnUlongIntDictionary(counts); } - ThreadStaticPools.ReturnUlongList(pending); - if (compositeStack is not null) + pending.Dispose(); + return result; + } + + private ref struct PendingBuffer + { + private Span _stack; + private ulong[]? _pooled; + private int _count; + + public PendingBuffer(Span stack) { - ThreadStaticPools.ReturnUlongStack(compositeStack); + _stack = stack; + _pooled = null; + _count = 0; } - return result; + public int Count => _count; + + public void Add(ulong value) + { + if (_pooled is null) + { + if (_count < _stack.Length) + { + _stack[_count] = value; + _count++; + return; + } + + int initialCapacity = _stack.Length == 0 ? 4 : _stack.Length * 2; + ulong[] newArray = ArrayPool.Shared.Rent(initialCapacity); + for (int i = 0; i < _count; i++) + { + newArray[i] = _stack[i]; + } + + _pooled = newArray; + } + + EnsurePooledCapacity(_count + 1); + _pooled![_count] = value; + _count++; + } + + private void EnsurePooledCapacity(int required) + { + ulong[] array = _pooled!; + if (required <= array.Length) + { + return; + } + + int newCapacity = array.Length * 2; + if (newCapacity < required) + { + newCapacity = required; + } + + ulong[] newArray = ArrayPool.Shared.Rent(newCapacity); + Array.Copy(array, 0, newArray, 0, _count); + ArrayPool.Shared.Return(array, clearArray: false); + _pooled = newArray; + } + + public void Clear() + { + _count = 0; + } + + public ReadOnlySpan AsSpan() + { + if (_pooled is null) + { + return _stack.Slice(0, _count); + } + + return new ReadOnlySpan(_pooled, 0, _count); + } + + public void Dispose() + { + if (_pooled is not null) + { + ArrayPool.Shared.Return(_pooled, clearArray: false); + _pooled = null; + } + + _count = 0; + } + } + + private ref struct CompositeStack + { + private Span _stack; + private ulong[]? _pooled; + private int _count; + + public CompositeStack(Span stack) + { + _stack = stack; + _pooled = null; + _count = 0; + } + + public int Count => _count; + + public void Push(ulong value) + { + if (_pooled is null) + { + if (_count < _stack.Length) + { + _stack[_count] = value; + _count++; + return; + } + + int initialCapacity = _stack.Length == 0 ? 4 : _stack.Length * 2; + ulong[] newArray = ArrayPool.Shared.Rent(initialCapacity); + for (int i = 0; i < _count; i++) + { + newArray[i] = _stack[i]; + } + + _pooled = newArray; + } + + EnsurePooledCapacity(_count + 1); + _pooled![_count] = value; + _count++; + } + + private void EnsurePooledCapacity(int required) + { + ulong[] array = _pooled!; + if (required <= array.Length) + { + return; + } + + int newCapacity = array.Length * 2; + if (newCapacity < required) + { + newCapacity = required; + } + + ulong[] newArray = ArrayPool.Shared.Rent(newCapacity); + Array.Copy(array, 0, newArray, 0, _count); + ArrayPool.Shared.Return(array, clearArray: false); + _pooled = newArray; + } + + public ulong Pop() + { + int index = _count - 1; + _count = index; + if (_pooled is null) + { + return _stack[index]; + } + + return _pooled[index]; + } + + public void Dispose() + { + if (_pooled is not null) + { + ArrayPool.Shared.Return(_pooled, clearArray: false); + _pooled = null; + } + + _count = 0; + } } private const int GpuSmallPrimeFactorSlots = 64; From 03751c27bedea82b3ee401b46b822280a9c84093 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 03:17:27 +0200 Subject: [PATCH 05/14] Add timed test mode for prime candidate scans --- EvenPerfectBitScanner/Program.cs | 143 ++++++++++++++++++++++++++++++- 1 file changed, 142 insertions(+), 1 deletion(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index b7e7ad9..b77ec55 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -1,5 +1,6 @@ using System.Buffers; using System.Collections.Generic; +using System.Diagnostics; using System.Globalization; using System.Numerics; using System.Runtime.CompilerServices; @@ -45,6 +46,11 @@ internal static class Program private static readonly Optimized PrimeIterator = new(); private static ulong _byDivisorStartPrime; + private static bool _testMode; + private static int _testTargetPrimeCount; + private static int _testProcessedPrimeCount; + private const string TestTimeFileName = "even_perfect_test_time.txt"; + [ThreadStatic] private static bool _lastCompositeP; @@ -95,6 +101,10 @@ internal static class Program ulong remainder = 0UL; bool startPrimeProvided = false; + _testMode = false; + _testTargetPrimeCount = 0; + _testProcessedPrimeCount = 0; + // NTT backend selection (GPU): reference vs staged for (; argIndex < args.Length; argIndex++) { @@ -201,6 +211,10 @@ internal static class Program { useModuloWorkaround = true; } + else if (arg.Equals("--test", StringComparison.OrdinalIgnoreCase)) + { + _testMode = true; + } else if (arg.Equals("--use-order", StringComparison.OrdinalIgnoreCase)) { useOrder = true; @@ -393,6 +407,20 @@ internal static class Program } } + if (_testMode) + { + currentP = 31UL; + remainder = currentP % 6UL; + useDivisor = false; + useByDivisor = false; + useResidue = false; + useLucas = false; + filterFile = string.Empty; + maxPrimeLimit = UInt128.MaxValue; + maxPrimeConfigured = false; + startPrimeProvided = true; + } + if (useByDivisor) { _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL; @@ -620,7 +648,7 @@ internal static class Program File.WriteAllText(ResultsFileName, $"p,searchedMersenne,detailedCheck,passedAllTests{Environment.NewLine}"); } - bool useFilter = !string.IsNullOrEmpty(filterFile); + bool useFilter = !_testMode && !string.IsNullOrEmpty(filterFile); HashSet filter = []; List byDivisorCandidates = []; ulong maxP = 0UL; @@ -737,7 +765,19 @@ internal static class Program } } + if (_testMode) + { + _testTargetPrimeCount = Math.Max(1, threadCount) * 3; + } + _state = ((long)currentP << 3) | (long)remainder; + + Stopwatch? stopwatch = null; + if (_testMode) + { + stopwatch = Stopwatch.StartNew(); + } + Task[] tasks = new Task[threadCount]; int taskIndex = 0; @@ -807,8 +847,23 @@ internal static class Program } Task.WaitAll(tasks); + + TimeSpan testElapsed = default; + bool reportTestTime = false; + if (stopwatch is not null) + { + stopwatch.Stop(); + testElapsed = stopwatch.Elapsed; + reportTestTime = true; + } + FlushBuffer(); StringBuilderPool.Return(_outputBuilder!); + + if (reportTestTime) + { + ReportTestTime(testElapsed); + } } private static List LoadByDivisorCandidates(string candidateFile, UInt128 maxPrimeLimit, bool maxPrimeConfigured, out int skippedByLimit) @@ -1053,6 +1108,21 @@ internal static class Program } } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void OnPrimeCandidateConfirmed() + { + if (!_testMode) + { + return; + } + + int processed = Interlocked.Increment(ref _testProcessedPrimeCount); + if (processed == _testTargetPrimeCount) + { + Volatile.Write(ref _limitReached, true); + } + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryAcquireConsoleSlot(bool primeCandidate, out bool primeFlag) { @@ -1189,6 +1259,75 @@ internal static class Program } } + private static void ReportTestTime(TimeSpan elapsed) + { + double elapsedSeconds = elapsed.TotalSeconds; + double elapsedMilliseconds = elapsed.TotalMilliseconds; + Console.WriteLine($"Test elapsed time: {elapsedSeconds:F3} s"); + + string filePath = GetTestTimeFilePath(); + bool hasPrevious = TryReadTestTime(filePath, out double previousMilliseconds); + string elapsedText = elapsedMilliseconds.ToString("F3", CultureInfo.InvariantCulture); + + if (!hasPrevious) + { + WriteTestTime(filePath, elapsedMilliseconds); + Console.WriteLine($"FIRST TIME: {elapsedText} ms"); + return; + } + + string previousText = previousMilliseconds.ToString("F3", CultureInfo.InvariantCulture); + if (elapsedMilliseconds < previousMilliseconds) + { + WriteTestTime(filePath, elapsedMilliseconds); + Console.WriteLine($"BETTER TIME: {elapsedText} ms (previous {previousText} ms)"); + } + else + { + Console.WriteLine($"WORSE TIME: {elapsedText} ms (best {previousText} ms)"); + } + } + + private static string GetTestTimeFilePath() + { + if (string.IsNullOrEmpty(_resultsDir)) + { + return TestTimeFileName; + } + + return Path.Combine(_resultsDir!, TestTimeFileName); + } + + private static bool TryReadTestTime(string filePath, out double milliseconds) + { + if (!File.Exists(filePath)) + { + milliseconds = 0.0; + return false; + } + + string content = File.ReadAllText(filePath).Trim(); + if (content.Length == 0) + { + milliseconds = 0.0; + return false; + } + + return double.TryParse(content, NumberStyles.Float, CultureInfo.InvariantCulture, out milliseconds); + } + + private static void WriteTestTime(string filePath, double milliseconds) + { + string? directory = Path.GetDirectoryName(filePath); + if (!string.IsNullOrEmpty(directory) && !Directory.Exists(directory)) + { + Directory.CreateDirectory(directory); + } + + string formatted = milliseconds.ToString("F6", CultureInfo.InvariantCulture); + File.WriteAllText(filePath, formatted); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static ulong CountOnes(ulong value) { @@ -1419,6 +1558,8 @@ internal static class Program _lastCompositeP = true; return false; } + + OnPrimeCandidateConfirmed(); } searchedMersenne = true; From ff724732f82b2e8cef0640f2065f478a99f09dba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 03:42:55 +0200 Subject: [PATCH 06/14] Allow test mode to seed by-divisor runs --- EvenPerfectBitScanner/Program.cs | 176 +++++++++++++++++++++---------- 1 file changed, 118 insertions(+), 58 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index b77ec55..1a31f3f 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -407,18 +407,17 @@ internal static class Program } } + int testPrimeCandidateLimit = 0; + if (_testMode) { currentP = 31UL; remainder = currentP % 6UL; - useDivisor = false; - useByDivisor = false; - useResidue = false; - useLucas = false; filterFile = string.Empty; maxPrimeLimit = UInt128.MaxValue; maxPrimeConfigured = false; startPrimeProvided = true; + testPrimeCandidateLimit = Math.Max(1, threadCount) * 3; } if (useByDivisor) @@ -430,7 +429,7 @@ internal static class Program _byDivisorStartPrime = 0UL; } - if (useByDivisor && string.IsNullOrEmpty(filterFile)) + if (useByDivisor && !_testMode && string.IsNullOrEmpty(filterFile)) { Console.WriteLine("--mersenne=bydivisor requires --filter-p=."); return; @@ -654,11 +653,16 @@ internal static class Program ulong maxP = 0UL; ulong[] localFilter = Array.Empty(); int filterCount = 0; - if (useFilter) + + if (_useByDivisorMode) { - Console.WriteLine("Loading filter..."); - if (_useByDivisorMode) + if (_testMode) { + byDivisorCandidates = BuildTestPrimeCandidates(testPrimeCandidateLimit); + } + else + { + Console.WriteLine("Loading filter..."); int skippedByLimit = 0; byDivisorCandidates = LoadByDivisorCandidates(filterFile, maxPrimeLimit, maxPrimeConfigured, out skippedByLimit); if (maxPrimeConfigured) @@ -675,63 +679,70 @@ internal static class Program } } } - else + } + else if (useFilter) + { + Console.WriteLine("Loading filter..."); + // TODO: Rent this filter buffer from ArrayPool and reuse it across reload batches so + // we do not allocate fresh arrays while replaying large result filters. + localFilter = new ulong[1024]; + filterCount = 0; + bool addedWithinLimit = false; + int skippedByLimit = 0; + LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) => { - // TODO: Rent this filter buffer from ArrayPool and reuse it across reload batches so - // we do not allocate fresh arrays while replaying large result filters. - localFilter = new ulong[1024]; - filterCount = 0; - bool addedWithinLimit = false; - int skippedByLimit = 0; - LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) => + if (!passedAllTests) { - if (!passedAllTests) - { - return; - } + return; + } - if (maxPrimeConfigured && (UInt128)p > maxPrimeLimit) - { - skippedByLimit++; - return; - } + if (maxPrimeConfigured && (UInt128)p > maxPrimeLimit) + { + skippedByLimit++; + return; + } - localFilter[filterCount++] = p; - addedWithinLimit = true; - if (p > maxP) - { - maxP = p; - } + localFilter[filterCount++] = p; + addedWithinLimit = true; + if (p > maxP) + { + maxP = p; + } - if (filterCount == 1024) - { - filter.AddRange(localFilter[..filterCount]); - filterCount = 0; - Console.WriteLine($"Added {p}"); - } - }); - - if (filterCount > 0) + if (filterCount == 1024) { filter.AddRange(localFilter[..filterCount]); + filterCount = 0; + Console.WriteLine($"Added {p}"); + } + }); + + if (filterCount > 0) + { + filter.AddRange(localFilter[..filterCount]); + } + + if (maxPrimeConfigured) + { + if (!addedWithinLimit) + { + Console.WriteLine("No filter candidates fall within the --max-prime limit."); + return; } - if (maxPrimeConfigured) + if (skippedByLimit > 0) { - if (!addedWithinLimit) - { - Console.WriteLine("No filter candidates fall within the --max-prime limit."); - return; - } - - if (skippedByLimit > 0) - { - Console.WriteLine($"Skipped {skippedByLimit} filter candidate(s) above --max-prime limit."); - } + Console.WriteLine($"Skipped {skippedByLimit} filter candidate(s) above --max-prime limit."); } } } + + if (_testMode) + { + _testTargetPrimeCount = testPrimeCandidateLimit; + } + if (_primeCount >= 2) { _primeFoundAfterInit = true; @@ -744,6 +755,12 @@ internal static class Program // Configure batch size for GPU primality sieve PrimeTester.GpuBatchSize = gpuPrimeBatch; + Stopwatch? stopwatch = null; + if (_testMode) + { + stopwatch = Stopwatch.StartNew(); + } + Console.WriteLine("Starting scan..."); if (_useByDivisorMode) @@ -751,6 +768,12 @@ internal static class Program RunByDivisorMode(byDivisorCandidates, threadCount); FlushBuffer(); StringBuilderPool.Return(_outputBuilder!); + if (stopwatch is not null) + { + stopwatch.Stop(); + ReportTestTime(stopwatch.Elapsed); + } + return; } @@ -767,17 +790,12 @@ internal static class Program if (_testMode) { - _testTargetPrimeCount = Math.Max(1, threadCount) * 3; + testPrimeCandidateLimit = Math.Max(1, threadCount) * 3; + _testTargetPrimeCount = testPrimeCandidateLimit; } _state = ((long)currentP << 3) | (long)remainder; - Stopwatch? stopwatch = null; - if (_testMode) - { - stopwatch = Stopwatch.StartNew(); - } - Task[] tasks = new Task[threadCount]; int taskIndex = 0; @@ -866,6 +884,48 @@ internal static class Program } } + private static List BuildTestPrimeCandidates(int targetCount) + { + if (targetCount <= 0) + { + return []; + } + + List candidates = new(targetCount); + ulong previous = 29UL; + int remaining = targetCount; + ulong nextPrime = 0UL; + + while (remaining > 0) + { + try + { + nextPrime = PrimeIterator.Next(in previous); + } + catch (InvalidOperationException) + { + break; + } + + if (nextPrime < 31UL) + { + previous = nextPrime; + continue; + } + + candidates.Add(nextPrime); + remaining--; + previous = nextPrime; + } + + if (remaining > 0) + { + Console.WriteLine("Unable to populate the requested number of test primes before reaching the 64-bit limit."); + } + + return candidates; + } + private static List LoadByDivisorCandidates(string candidateFile, UInt128 maxPrimeLimit, bool maxPrimeConfigured, out int skippedByLimit) { List candidates = []; From 3e7ee6ab464790f9cc0796bfcb4518c09aee9760 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 04:52:17 +0200 Subject: [PATCH 07/14] Restore pooled buffers in partial factorization --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 197 ++------------------ 1 file changed, 13 insertions(+), 184 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index da15d5c..b3da51d 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -799,8 +799,8 @@ internal static partial class PrimeOrderCalculator gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); } - Span pendingStorage = stackalloc ulong[4]; - PendingBuffer pending = new PendingBuffer(pendingStorage); + List pending = ThreadStaticPools.RentUlongList(2); + Stack? compositeStack = null; PartialFactorResult result; if (!gpuFactored) @@ -822,10 +822,8 @@ internal static partial class PrimeOrderCalculator if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) { - Span compositeStorage = stackalloc ulong[8]; - CompositeStack compositeStack = new CompositeStack(compositeStorage); - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count * 2, 4)); compositeStack.Push(remaining); pending.Clear(); @@ -862,16 +860,13 @@ internal static partial class PrimeOrderCalculator compositeStack.Push(factor); compositeStack.Push(quotient); } - - compositeStack.Dispose(); } ulong cofactor = 1UL; - ReadOnlySpan pendingSpan = pending.AsSpan(); - int pendingCount = pendingSpan.Length; + int pendingCount = pending.Count; for (int i = 0; i < pendingCount; i++) { - ulong composite = pendingSpan[i]; + ulong composite = pending[i]; bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) { @@ -961,183 +956,17 @@ internal static partial class PrimeOrderCalculator ThreadStaticPools.ReturnUlongIntDictionary(counts); } - pending.Dispose(); + if (compositeStack is not null) + { + compositeStack.Clear(); + ThreadStaticPools.ReturnUlongStack(compositeStack); + } + + pending.Clear(); + ThreadStaticPools.ReturnUlongList(pending); return result; } - private ref struct PendingBuffer - { - private Span _stack; - private ulong[]? _pooled; - private int _count; - - public PendingBuffer(Span stack) - { - _stack = stack; - _pooled = null; - _count = 0; - } - - public int Count => _count; - - public void Add(ulong value) - { - if (_pooled is null) - { - if (_count < _stack.Length) - { - _stack[_count] = value; - _count++; - return; - } - - int initialCapacity = _stack.Length == 0 ? 4 : _stack.Length * 2; - ulong[] newArray = ArrayPool.Shared.Rent(initialCapacity); - for (int i = 0; i < _count; i++) - { - newArray[i] = _stack[i]; - } - - _pooled = newArray; - } - - EnsurePooledCapacity(_count + 1); - _pooled![_count] = value; - _count++; - } - - private void EnsurePooledCapacity(int required) - { - ulong[] array = _pooled!; - if (required <= array.Length) - { - return; - } - - int newCapacity = array.Length * 2; - if (newCapacity < required) - { - newCapacity = required; - } - - ulong[] newArray = ArrayPool.Shared.Rent(newCapacity); - Array.Copy(array, 0, newArray, 0, _count); - ArrayPool.Shared.Return(array, clearArray: false); - _pooled = newArray; - } - - public void Clear() - { - _count = 0; - } - - public ReadOnlySpan AsSpan() - { - if (_pooled is null) - { - return _stack.Slice(0, _count); - } - - return new ReadOnlySpan(_pooled, 0, _count); - } - - public void Dispose() - { - if (_pooled is not null) - { - ArrayPool.Shared.Return(_pooled, clearArray: false); - _pooled = null; - } - - _count = 0; - } - } - - private ref struct CompositeStack - { - private Span _stack; - private ulong[]? _pooled; - private int _count; - - public CompositeStack(Span stack) - { - _stack = stack; - _pooled = null; - _count = 0; - } - - public int Count => _count; - - public void Push(ulong value) - { - if (_pooled is null) - { - if (_count < _stack.Length) - { - _stack[_count] = value; - _count++; - return; - } - - int initialCapacity = _stack.Length == 0 ? 4 : _stack.Length * 2; - ulong[] newArray = ArrayPool.Shared.Rent(initialCapacity); - for (int i = 0; i < _count; i++) - { - newArray[i] = _stack[i]; - } - - _pooled = newArray; - } - - EnsurePooledCapacity(_count + 1); - _pooled![_count] = value; - _count++; - } - - private void EnsurePooledCapacity(int required) - { - ulong[] array = _pooled!; - if (required <= array.Length) - { - return; - } - - int newCapacity = array.Length * 2; - if (newCapacity < required) - { - newCapacity = required; - } - - ulong[] newArray = ArrayPool.Shared.Rent(newCapacity); - Array.Copy(array, 0, newArray, 0, _count); - ArrayPool.Shared.Return(array, clearArray: false); - _pooled = newArray; - } - - public ulong Pop() - { - int index = _count - 1; - _count = index; - if (_pooled is null) - { - return _stack[index]; - } - - return _pooled[index]; - } - - public void Dispose() - { - if (_pooled is not null) - { - ArrayPool.Shared.Return(_pooled, clearArray: false); - _pooled = null; - } - - _count = 0; - } - } - private const int GpuSmallPrimeFactorSlots = 64; private static bool TryPopulateSmallPrimeFactorsCpu( From 949b5f891c6b39d3ac3310e970a67ba84f8468e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 04:52:22 +0200 Subject: [PATCH 08/14] Return windowed Montgomery odd powers to pooling --- PerfectNumbers.Core/ULongExtensions.cs | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 0271021..c09ce8e 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -346,7 +346,8 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - Span oddPowers = stackalloc ulong[oddPowerCount]; + ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers); int index = bitLength - 1; @@ -385,13 +386,14 @@ public static class ULongExtensions index = windowStart - 1; } + ThreadStaticPools.UlongPool.Return(oddPowersArray); - if (keepMontgomery) + if (!keepMontgomery) { - return result; + result = result.MontgomeryMultiply(1UL, modulus, nPrime); } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result; } [MethodImpl(MethodImplOptions.AggressiveInlining)] From fb08a61c66c81eeacfb5ceb021092e2860fdeea4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 17:18:19 +0200 Subject: [PATCH 09/14] Restore pooled small factor buffers --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 63 ++++++++++++--------- PerfectNumbers.Core/ThreadStaticPools.cs | 12 ++++ 2 files changed, 49 insertions(+), 26 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index b3da51d..5ca1876 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -781,8 +781,10 @@ internal static partial class PrimeOrderCalculator } const int FactorSlotCount = GpuSmallPrimeFactorSlots; - Span primeSlots = stackalloc ulong[FactorSlotCount]; - Span exponentSlots = stackalloc int[FactorSlotCount]; + ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); + int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); + Span primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount); + Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); primeSlots.Clear(); exponentSlots.Clear(); @@ -951,6 +953,9 @@ internal static partial class PrimeOrderCalculator result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); ReturnResult: + ThreadStaticPools.UlongPool.Return(primeSlotsArray); + ThreadStaticPools.IntPool.Return(exponentSlotsArray); + if (counts is not null) { ThreadStaticPools.ReturnUlongIntDictionary(counts); @@ -1079,34 +1084,40 @@ internal static partial class PrimeOrderCalculator 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(); + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) + { + ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); + int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); + Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + primeBuffer.Clear(); + exponentBuffer.Clear(); - if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) - { - remaining = value; - return false; - } + if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) + { + remaining = value; + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + 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; - } + 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; - } + counts[primeValue] = exponent; + } - return true; - } + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + return true; + } private static long CreateDeadlineTimestamp(int milliseconds) { if (milliseconds <= 0) diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 9ba6060..ed00820 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -79,6 +79,18 @@ namespace PerfectNumbers.Core } } + [ThreadStatic] + private static ArrayPool? _intPool; + + public static ArrayPool IntPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _intPool ??= ArrayPool.Create(); + } + } + [ThreadStatic] private static List>? _ulongListPool; From e53bb6ba43772fad941c176b748c447fb38673f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 16 Oct 2025 17:18:27 +0200 Subject: [PATCH 10/14] Track small-prime remainders with cycle steppers --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 109 +++++++++++++++++++- 1 file changed, 108 insertions(+), 1 deletion(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 5ca1876..0aaca1d 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -47,6 +47,18 @@ internal static partial class PrimeOrderCalculator [ThreadStatic] private static bool s_debugLoggingEnabled; + [ThreadStatic] + private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers; + + [ThreadStatic] + private static ulong[]? s_smallPrimeStepperPreviousValues; + + [ThreadStatic] + private static ulong[]? s_smallPrimeStepperRemainders; + + [ThreadStatic] + private static bool[]? s_smallPrimeStepperHasState; + @@ -972,6 +984,62 @@ internal static partial class PrimeOrderCalculator return result; } + private static CycleRemainderStepper[] EnsureSmallPrimeCycleSteppers(uint[] primes) + { + int primeCount = primes.Length; + CycleRemainderStepper[]? steppers = s_smallPrimeCycleSteppers; + if (steppers is null || steppers.Length < primeCount) + { + int length = primeCount; + steppers = new CycleRemainderStepper[length]; + for (int i = 0; i < length; i++) + { + steppers[i] = new CycleRemainderStepper(primes[i]); + } + + s_smallPrimeCycleSteppers = steppers; + s_smallPrimeStepperPreviousValues = new ulong[length]; + s_smallPrimeStepperRemainders = new ulong[length]; + s_smallPrimeStepperHasState = new bool[length]; + return steppers; + } + + if (s_smallPrimeStepperPreviousValues is null || s_smallPrimeStepperPreviousValues.Length != steppers.Length) + { + var previousValues = new ulong[steppers.Length]; + if (s_smallPrimeStepperPreviousValues is { } existingPrevious) + { + System.Array.Copy(existingPrevious, previousValues, Math.Min(existingPrevious.Length, previousValues.Length)); + } + + s_smallPrimeStepperPreviousValues = previousValues; + } + + if (s_smallPrimeStepperRemainders is null || s_smallPrimeStepperRemainders.Length != steppers.Length) + { + var remainders = new ulong[steppers.Length]; + if (s_smallPrimeStepperRemainders is { } existingRemainders) + { + System.Array.Copy(existingRemainders, remainders, Math.Min(existingRemainders.Length, remainders.Length)); + } + + s_smallPrimeStepperRemainders = remainders; + } + + if (s_smallPrimeStepperHasState is null || s_smallPrimeStepperHasState.Length != steppers.Length) + { + var hasState = new bool[steppers.Length]; + if (s_smallPrimeStepperHasState is { } existingHasState) + { + System.Array.Copy(existingHasState, hasState, Math.Min(existingHasState.Length, hasState.Length)); + } + + s_smallPrimeStepperHasState = hasState; + } + + return steppers; + } + private const int GpuSmallPrimeFactorSlots = 64; private static bool TryPopulateSmallPrimeFactorsCpu( @@ -1001,6 +1069,11 @@ internal static partial class PrimeOrderCalculator int primeCount = primes.Length; ulong remainingLocal = value; uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + CycleRemainderStepper[] cycleSteppers = EnsureSmallPrimeCycleSteppers(primes); + ulong[] stepperPreviousValues = s_smallPrimeStepperPreviousValues!; + ulong[] stepperRemainders = s_smallPrimeStepperRemainders!; + bool[] stepperHasState = s_smallPrimeStepperHasState!; + ulong valueForSteppers = value; for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) { @@ -1017,11 +1090,45 @@ internal static partial class PrimeOrderCalculator } ulong primeValue = primeCandidate; - if ((remainingLocal % primeValue) != 0UL) + if (primeValue == 0UL) { continue; } + ref CycleRemainderStepper stepper = ref cycleSteppers[i]; + ulong remainder; + + if (!stepperHasState[i]) + { + stepper.Reset(); + remainder = stepper.Initialize(valueForSteppers); + stepperHasState[i] = true; + } + else + { + ulong previousValue = stepperPreviousValues[i]; + if (valueForSteppers > previousValue) + { + remainder = stepper.ComputeNext(valueForSteppers); + } + else if (valueForSteppers < previousValue) + { + stepper.Reset(); + remainder = stepper.Initialize(valueForSteppers); + } + else + { + remainder = stepperRemainders[i]; + } + } + + stepperPreviousValues[i] = valueForSteppers; + stepperRemainders[i] = remainder; + + if (remainder != 0UL) + { + continue; + } int exponent = 0; do { From 07395830a5a31e7300933c1507fb8d713a99d991 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 16 Oct 2025 18:13:28 +0200 Subject: [PATCH 11/14] Unnecessary cleanup commented out --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 0aaca1d..0829ff9 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -797,8 +797,9 @@ internal static partial class PrimeOrderCalculator int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); Span primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount); Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); - primeSlots.Clear(); - exponentSlots.Clear(); + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); int factorCount = 0; Dictionary? counts = null; @@ -819,8 +820,9 @@ internal static partial class PrimeOrderCalculator if (!gpuFactored) { - primeSlots.Clear(); - exponentSlots.Clear(); + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) { useDictionary = true; From aa4ff799f34bd5883aa97bef28533419a33d327e Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 16 Oct 2025 19:46:25 +0200 Subject: [PATCH 12/14] Small optimizations --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 843 ++++++++++---------- 1 file changed, 428 insertions(+), 415 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 0829ff9..89efdc1 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -47,17 +47,17 @@ internal static partial class PrimeOrderCalculator [ThreadStatic] private static bool s_debugLoggingEnabled; - [ThreadStatic] - private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers; + [ThreadStatic] + private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers; - [ThreadStatic] - private static ulong[]? s_smallPrimeStepperPreviousValues; + [ThreadStatic] + private static ulong[]? s_smallPrimeStepperPreviousValues; - [ThreadStatic] - private static ulong[]? s_smallPrimeStepperRemainders; + [ThreadStatic] + private static ulong[]? s_smallPrimeStepperRemainders; - [ThreadStatic] - private static bool[]? s_smallPrimeStepperHasState; + [ThreadStatic] + private static bool[]? s_smallPrimeStepperHasState; @@ -531,7 +531,7 @@ internal static partial class PrimeOrderCalculator } ThreadStaticPools.ReturnUlongList(candidates); - // DebugLog("No candidate confirmed"); + // DebugLog("No candidate confirmed"); return false; } finally @@ -786,374 +786,387 @@ internal static partial class PrimeOrderCalculator } private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } - const int FactorSlotCount = GpuSmallPrimeFactorSlots; - ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); - int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); - Span primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount); - Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); + const int FactorSlotCount = GpuSmallPrimeFactorSlots; - int factorCount = 0; - Dictionary? counts = null; - bool useDictionary = false; + // stackalloc is faster than pooling + // ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); + // int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); + // Span primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount); + // Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); + Span primeSlots = stackalloc ulong[FactorSlotCount]; + Span exponentSlots = stackalloc int[FactorSlotCount]; + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; + int factorCount = 0; + Dictionary? counts = null; + bool useDictionary = false; - bool gpuFactored = false; - if (IsGpuHeuristicDevice) - { - gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); - } + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong remaining = value; - List pending = ThreadStaticPools.RentUlongList(2); - Stack? compositeStack = null; - PartialFactorResult result; + bool gpuFactored = false; + if (IsGpuHeuristicDevice) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); + } - if (!gpuFactored) - { - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); - if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) - { - useDictionary = true; - counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(factorCount, 8)); - remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); - } - } + List pending = ThreadStaticPools.RentUlongList(2); + Stack? compositeStack = null; + PartialFactorResult result; - if (remaining > 1UL) - { - pending.Add(remaining); - } + if (!gpuFactored) + { + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); + if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) + { + useDictionary = true; + counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(factorCount, 8)); + remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + } + } - if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count * 2, 4)); - compositeStack.Push(remaining); - pending.Clear(); + if (remaining > 1UL) + { + pending.Add(remaining); + } - long timestamp = 0L; - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - if (composite <= 1UL) - { - continue; - } + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count << 1, 4)); + compositeStack.Push(remaining); + pending.Clear(); - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - if (isPrime) - { - AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - continue; - } + long timestamp = 0L; + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + // composite will never be smaller than 1 on the execution path + // if (composite <= 1UL) + // { + // continue; + // } - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - pending.Add(composite); - 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 (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(composite); - continue; - } + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + pending.Add(composite); + continue; + } - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - } + if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(composite); + continue; + } - 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); - } - } + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + } - ArrayPool pool = ThreadStaticPools.FactorEntryPool; - if (useDictionary) - { - if ((counts is null || counts.Count == 0) && cofactor == value) - { - result = PartialFactorResult.Rent(null, value, false, 0); - goto ReturnResult; - } + 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); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * composite); + } + } - if (counts is null || counts.Count == 0) - { - result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - goto ReturnResult; - } + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + if (useDictionary) + { + if ((counts is null || counts.Count == 0) && cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; + } - FactorEntry[] factors = pool.Rent(counts.Count); - int index = 0; - foreach (KeyValuePair entry in counts) - { - factors[index] = new FactorEntry(entry.Key, entry.Value); - index++; - } + if (counts is null || counts.Count == 0) + { + result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + goto ReturnResult; + } - Span span = factors.AsSpan(0, index); - span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + FactorEntry[] factors = pool.Rent(counts.Count); + int index = 0; + foreach (KeyValuePair entry in counts) + { + factors[index] = new FactorEntry(entry.Key, entry.Value); + index++; + } - bool fullyFactored = cofactor == 1UL; - result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); - goto ReturnResult; - } + Span span = factors.AsSpan(0, index); + span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - int actualCount = 0; - for (int i = 0; i < factorCount; i++) - { - if (primeSlots[i] != 0UL && exponentSlots[i] != 0) - { - actualCount++; - } - } + bool fullyFactored = cofactor == 1UL; + result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); + goto ReturnResult; + } - if (actualCount == 0) - { - if (cofactor == value) - { - result = PartialFactorResult.Rent(null, value, false, 0); - goto ReturnResult; - } + int actualCount = 0; + for (int i = 0; i < factorCount; i++) + { + if (primeSlots[i] != 0UL && exponentSlots[i] != 0) + { + actualCount++; + } + } - result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - goto ReturnResult; - } + if (actualCount == 0) + { + if (cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; + } - 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; - } + result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + goto ReturnResult; + } - array[arrayIndex] = new FactorEntry(primeValue, exponentValue); - arrayIndex++; - } + 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; + } - Span arraySpan = array.AsSpan(0, arrayIndex); - arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - bool fullyFactoredArray = cofactor == 1UL; - result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); + array[arrayIndex] = new FactorEntry(primeValue, exponentValue); + arrayIndex++; + } - ReturnResult: - ThreadStaticPools.UlongPool.Return(primeSlotsArray); - ThreadStaticPools.IntPool.Return(exponentSlotsArray); + Span arraySpan = array.AsSpan(0, arrayIndex); + arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + bool fullyFactoredArray = cofactor == 1UL; + result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); - if (counts is not null) - { - ThreadStaticPools.ReturnUlongIntDictionary(counts); - } + ReturnResult: + // ThreadStaticPools.UlongPool.Return(primeSlotsArray); + // ThreadStaticPools.IntPool.Return(exponentSlotsArray); - if (compositeStack is not null) - { - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); - } + if (counts is not null) + { + ThreadStaticPools.ReturnUlongIntDictionary(counts); + } - pending.Clear(); - ThreadStaticPools.ReturnUlongList(pending); - return result; - } + if (compositeStack is not null) + { + compositeStack.Clear(); + ThreadStaticPools.ReturnUlongStack(compositeStack); + } - private static CycleRemainderStepper[] EnsureSmallPrimeCycleSteppers(uint[] primes) - { - int primeCount = primes.Length; - CycleRemainderStepper[]? steppers = s_smallPrimeCycleSteppers; - if (steppers is null || steppers.Length < primeCount) - { - int length = primeCount; - steppers = new CycleRemainderStepper[length]; - for (int i = 0; i < length; i++) - { - steppers[i] = new CycleRemainderStepper(primes[i]); - } + pending.Clear(); + ThreadStaticPools.ReturnUlongList(pending); + return result; + } - s_smallPrimeCycleSteppers = steppers; - s_smallPrimeStepperPreviousValues = new ulong[length]; - s_smallPrimeStepperRemainders = new ulong[length]; - s_smallPrimeStepperHasState = new bool[length]; - return steppers; - } + private static CycleRemainderStepper[] EnsureSmallPrimeCycleSteppers(uint[] primes) + { + int primeCount = primes.Length; + CycleRemainderStepper[]? steppers = s_smallPrimeCycleSteppers; + if (steppers is null || steppers.Length < primeCount) + { + int length = primeCount; + steppers = new CycleRemainderStepper[length]; + for (int i = 0; i < length; i++) + { + steppers[i] = new CycleRemainderStepper(primes[i]); + } - if (s_smallPrimeStepperPreviousValues is null || s_smallPrimeStepperPreviousValues.Length != steppers.Length) - { - var previousValues = new ulong[steppers.Length]; - if (s_smallPrimeStepperPreviousValues is { } existingPrevious) - { - System.Array.Copy(existingPrevious, previousValues, Math.Min(existingPrevious.Length, previousValues.Length)); - } + s_smallPrimeCycleSteppers = steppers; + s_smallPrimeStepperPreviousValues = new ulong[length]; + s_smallPrimeStepperRemainders = new ulong[length]; + s_smallPrimeStepperHasState = new bool[length]; + return steppers; + } - s_smallPrimeStepperPreviousValues = previousValues; - } + if (s_smallPrimeStepperPreviousValues is null || s_smallPrimeStepperPreviousValues.Length != steppers.Length) + { + var previousValues = new ulong[steppers.Length]; + if (s_smallPrimeStepperPreviousValues is { } existingPrevious) + { + System.Array.Copy(existingPrevious, previousValues, Math.Min(existingPrevious.Length, previousValues.Length)); + } - if (s_smallPrimeStepperRemainders is null || s_smallPrimeStepperRemainders.Length != steppers.Length) - { - var remainders = new ulong[steppers.Length]; - if (s_smallPrimeStepperRemainders is { } existingRemainders) - { - System.Array.Copy(existingRemainders, remainders, Math.Min(existingRemainders.Length, remainders.Length)); - } + s_smallPrimeStepperPreviousValues = previousValues; + } - s_smallPrimeStepperRemainders = remainders; - } + if (s_smallPrimeStepperRemainders is null || s_smallPrimeStepperRemainders.Length != steppers.Length) + { + var remainders = new ulong[steppers.Length]; + if (s_smallPrimeStepperRemainders is { } existingRemainders) + { + System.Array.Copy(existingRemainders, remainders, Math.Min(existingRemainders.Length, remainders.Length)); + } - if (s_smallPrimeStepperHasState is null || s_smallPrimeStepperHasState.Length != steppers.Length) - { - var hasState = new bool[steppers.Length]; - if (s_smallPrimeStepperHasState is { } existingHasState) - { - System.Array.Copy(existingHasState, hasState, Math.Min(existingHasState.Length, hasState.Length)); - } + s_smallPrimeStepperRemainders = remainders; + } - s_smallPrimeStepperHasState = hasState; - } + if (s_smallPrimeStepperHasState is null || s_smallPrimeStepperHasState.Length != steppers.Length) + { + var hasState = new bool[steppers.Length]; + if (s_smallPrimeStepperHasState is { } existingHasState) + { + System.Array.Copy(existingHasState, hasState, Math.Min(existingHasState.Length, hasState.Length)); + } - return steppers; - } + s_smallPrimeStepperHasState = hasState; + } - private const int GpuSmallPrimeFactorSlots = 64; + return steppers; + } - private static bool TryPopulateSmallPrimeFactorsCpu( - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining) - { - factorCount = 0; - remaining = value; + private const int GpuSmallPrimeFactorSlots = 64; - if (value <= 1UL) - { - return true; - } + private static bool TryPopulateSmallPrimeFactorsCpu( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) + { + factorCount = 0; + remaining = value; - int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); - if (capacity == 0) - { - return false; - } + // value will never <= 1 in production code + // if (value <= 1UL) + // { + // return true; + // } - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; - ulong remainingLocal = value; - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - CycleRemainderStepper[] cycleSteppers = EnsureSmallPrimeCycleSteppers(primes); - ulong[] stepperPreviousValues = s_smallPrimeStepperPreviousValues!; - ulong[] stepperRemainders = s_smallPrimeStepperRemainders!; - bool[] stepperHasState = s_smallPrimeStepperHasState!; - ulong valueForSteppers = value; + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + // capacity will never equal 0 in production code + // if (capacity == 0) + // { + // return false; + // } - for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > effectiveLimit) - { - break; - } + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + ulong remainingLocal = value; + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + CycleRemainderStepper[] cycleSteppers = EnsureSmallPrimeCycleSteppers(primes); + ulong[] stepperPreviousValues = s_smallPrimeStepperPreviousValues!; + ulong[] stepperRemainders = s_smallPrimeStepperRemainders!; + bool[] stepperHasState = s_smallPrimeStepperHasState!; + ulong valueForSteppers = value; - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remainingLocal) - { - break; - } + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } - ulong primeValue = primeCandidate; - if (primeValue == 0UL) - { - continue; - } + ulong primeSquare = squares[i]; + // primeSquare will never == 0 in production code + // if (primeSquare != 0 && primeSquare > remainingLocal) + if (primeSquare > remainingLocal) + { + break; + } - ref CycleRemainderStepper stepper = ref cycleSteppers[i]; - ulong remainder; + // primeCandidate will never equal 0 in production code + // if (primeCandidate == 0UL) + // { + // continue; + // } - if (!stepperHasState[i]) - { - stepper.Reset(); - remainder = stepper.Initialize(valueForSteppers); - stepperHasState[i] = true; - } - else - { - ulong previousValue = stepperPreviousValues[i]; - if (valueForSteppers > previousValue) - { - remainder = stepper.ComputeNext(valueForSteppers); - } - else if (valueForSteppers < previousValue) - { - stepper.Reset(); - remainder = stepper.Initialize(valueForSteppers); - } - else - { - remainder = stepperRemainders[i]; - } - } + ref CycleRemainderStepper stepper = ref cycleSteppers[i]; + ulong remainder; - stepperPreviousValues[i] = valueForSteppers; - stepperRemainders[i] = remainder; + if (!stepperHasState[i]) + { + stepper.Reset(); + remainder = stepper.Initialize(valueForSteppers); + stepperHasState[i] = true; + } + else + { + ulong previousValue = stepperPreviousValues[i]; + if (valueForSteppers > previousValue) + { + remainder = stepper.ComputeNext(valueForSteppers); + } + else if (valueForSteppers < previousValue) + { + stepper.Reset(); + remainder = stepper.Initialize(valueForSteppers); + } + else + { + remainder = stepperRemainders[i]; + } + } - if (remainder != 0UL) - { - continue; - } - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); + stepperPreviousValues[i] = valueForSteppers; + stepperRemainders[i] = remainder; - if (factorCount >= capacity) - { - return false; - } + if (remainder != 0UL) + { + continue; + } - primeTargets[factorCount] = primeValue; - exponentTargets[factorCount] = exponent; - factorCount++; - } + // TODO: Can we calculate this using residue steppers? We can try reusing modulo calculation for division result or the other way around + int exponent = 0; + do + { + remainingLocal /= primeCandidate; + exponent++; + } + while ((remainingLocal % primeCandidate) == 0UL); - remaining = remainingLocal; - return true; - } + if (factorCount >= capacity) + { + throw new InvalidOperationException($"Capacity is smaller than factor count"); + } - private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) + primeTargets[factorCount] = primeCandidate; + 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; @@ -1193,111 +1206,111 @@ internal static partial class PrimeOrderCalculator return remaining; } - private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) - { - ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); - int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); - Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - primeBuffer.Clear(); - exponentBuffer.Clear(); + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) + { + ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); + int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); + Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + primeBuffer.Clear(); + exponentBuffer.Clear(); - if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) - { - remaining = value; - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - return false; - } + if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) + { + remaining = value; + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + 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; - } + 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; - } + counts[primeValue] = exponent; + } - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - return true; - } - private static long CreateDeadlineTimestamp(int milliseconds) - { - if (milliseconds <= 0) - { - return long.MaxValue; - } + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + return true; + } + private static long CreateDeadlineTimestamp(int milliseconds) + { + if (milliseconds <= 0) + { + return long.MaxValue; + } - long stopwatchTicks = ConvertMillisecondsToStopwatchTicks(milliseconds); - long startTimestamp = Stopwatch.GetTimestamp(); - long deadline = startTimestamp + stopwatchTicks; - if (deadline < startTimestamp) - { - return long.MaxValue; - } + long stopwatchTicks = ConvertMillisecondsToStopwatchTicks(milliseconds); + long startTimestamp = Stopwatch.GetTimestamp(); + long deadline = startTimestamp + stopwatchTicks; + if (deadline < startTimestamp) + { + return long.MaxValue; + } - return deadline; - } + return deadline; + } - private static long ConvertMillisecondsToStopwatchTicks(int milliseconds) - { - long frequency = Stopwatch.Frequency; - long baseTicks = (frequency / 1000L) * milliseconds; - long remainder = (frequency % 1000L) * milliseconds; - return baseTicks + (remainder / 1000L); - } + private static long ConvertMillisecondsToStopwatchTicks(int milliseconds) + { + long frequency = Stopwatch.Frequency; + long baseTicks = (frequency / 1000L) * milliseconds; + long remainder = (frequency % 1000L) * milliseconds; + return baseTicks + (remainder / 1000L); + } - private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor) - { - factor = 0UL; - if ((n & 1UL) == 0UL) - { - factor = 2UL; - return true; - } + private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor) + { + factor = 0UL; + if ((n & 1UL) == 0UL) + { + factor = 2UL; + return true; + } - long timestamp = 0L; // reused for deadline checks. - while (true) - { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } + long timestamp = 0L; // reused for deadline checks. + while (true) + { + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + return false; + } - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; + ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; - while (d == 1UL) - { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } + while (d == 1UL) + { + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + 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); - } + 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; - } - } - } + if (d != n) + { + factor = d; + return true; + } + } + } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) From 10751bf8418872b570ec8f02da76dcbc684f8b13 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 16 Oct 2025 21:10:58 +0200 Subject: [PATCH 13/14] Concurrent dictionary added for IsPrime results caching to utilize the repetition of composites. --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 22 +++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 89efdc1..fc4b039 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -6,6 +6,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using PerfectNumbers.Core.Gpu; +using System.Collections.Concurrent; namespace PerfectNumbers.Core; @@ -785,8 +786,15 @@ internal static partial class PrimeOrderCalculator return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; } + // [ThreadStatic] + private static ConcurrentDictionary? _primalityCache = []; + + // [ThreadStatic] + // private static UInt128 _isPrimeHits; + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) { + var primalityCache = _primalityCache;// ??= []; if (value <= 1UL) { return PartialFactorResult.Empty; @@ -857,7 +865,17 @@ internal static partial class PrimeOrderCalculator // continue; // } - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (!primalityCache!.TryGetValue(composite, out bool isPrime)) + { + // Console.WriteLine($"Checking composite {composite}"); + isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + _ = primalityCache.TryAdd(composite, isPrime); + } + // else + // { + // Console.WriteLine($"isPrimeHits = {++_isPrimeHits}"); + // } + if (isPrime) { AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); @@ -1143,7 +1161,7 @@ internal static partial class PrimeOrderCalculator continue; } - // TODO: Can we calculate this using residue steppers? We can try reusing modulo calculation for division result or the other way around + // TODO: Can we calculate this using residue stepping? We can try reusing modulo calculation for division result or the other way around int exponent = 0; do { From 4c81f9e5a2e28eddd632831c1dc502eb48fce617 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 16 Oct 2025 21:21:40 +0200 Subject: [PATCH 14/14] IsPrime cache reused --- PerfectNumbers.Core/PrimeOrderCalculator.cs | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index fc4b039..68fe8a9 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -906,7 +906,12 @@ internal static partial class PrimeOrderCalculator for (int i = 0; i < pendingCount; i++) { ulong composite = pending[i]; - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (!primalityCache!.TryGetValue(composite, out bool isPrime)) + { + // Console.WriteLine($"Checking composite {composite}"); + isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + _ = primalityCache.TryAdd(composite, isPrime); + } if (isPrime) { AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);