diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index b7e7ad9..1a31f3f 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,19 @@ internal static class Program } } + int testPrimeCandidateLimit = 0; + + if (_testMode) + { + currentP = 31UL; + remainder = currentP % 6UL; + filterFile = string.Empty; + maxPrimeLimit = UInt128.MaxValue; + maxPrimeConfigured = false; + startPrimeProvided = true; + testPrimeCandidateLimit = Math.Max(1, threadCount) * 3; + } + if (useByDivisor) { _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL; @@ -402,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; @@ -620,17 +647,22 @@ 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; 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) @@ -647,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; @@ -716,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) @@ -723,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; } @@ -737,7 +788,14 @@ internal static class Program } } + if (_testMode) + { + testPrimeCandidateLimit = Math.Max(1, threadCount) * 3; + _testTargetPrimeCount = testPrimeCandidateLimit; + } + _state = ((long)currentP << 3) | (long)remainder; + Task[] tasks = new Task[threadCount]; int taskIndex = 0; @@ -807,8 +865,65 @@ 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 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) @@ -1053,6 +1168,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 +1319,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 +1618,8 @@ internal static class Program _lastCompositeP = true; return false; } + + OnPrimeCandidateConfirmed(); } searchedMersenne = true; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 39d710c..68fe8a9 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -1,10 +1,12 @@ using System.Buffers; +using System.Collections.Generic; using System.Buffers.Binary; using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using PerfectNumbers.Core.Gpu; +using System.Collections.Concurrent; namespace PerfectNumbers.Core; @@ -46,15 +48,22 @@ internal static partial class PrimeOrderCalculator [ThreadStatic] private static bool s_debugLoggingEnabled; + [ThreadStatic] + private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers; [ThreadStatic] - private static List? s_partialFactorPending; + private static ulong[]? s_smallPrimeStepperPreviousValues; [ThreadStatic] - private static Stack? s_partialFactorCompositeStack; + private static ulong[]? s_smallPrimeStepperRemainders; [ThreadStatic] - private static Dictionary? s_partialFactorCounts; + private static bool[]? s_smallPrimeStepperHasState; + + + + + private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; @@ -402,12 +411,13 @@ internal static partial class PrimeOrderCalculator } int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; - List candidates = new(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) { + ThreadStaticPools.ReturnUlongList(candidates); return false; } @@ -508,6 +518,7 @@ internal static partial class PrimeOrderCalculator pool.Return(gpuPool, clearArray: false); } + ThreadStaticPools.ReturnUlongList(candidates); result = candidate; return true; } @@ -520,6 +531,7 @@ internal static partial class PrimeOrderCalculator index += batchSize; } + ThreadStaticPools.ReturnUlongList(candidates); // DebugLog("No candidate confirmed"); return false; } @@ -774,67 +786,32 @@ 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; - } + // [ThreadStatic] + private static ConcurrentDictionary? _primalityCache = []; - 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; - } + // [ThreadStatic] + // private static UInt128 _isPrimeHits; private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) { + var primalityCache = _primalityCache;// ??= []; if (value <= 1UL) { return PartialFactorResult.Empty; } const int FactorSlotCount = GpuSmallPrimeFactorSlots; + + // 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]; - 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; @@ -849,190 +826,288 @@ internal static partial class PrimeOrderCalculator gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); } - List pending = AcquirePendingCompositeList(2); + List pending = ThreadStaticPools.RentUlongList(2); Stack? compositeStack = null; + PartialFactorResult result; - try + if (!gpuFactored) { - 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)) { - 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 (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(); + + 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; + // } + + if (!primalityCache!.TryGetValue(composite, out bool isPrime)) { - useDictionary = true; - counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8)); - remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + // Console.WriteLine($"Checking composite {composite}"); + isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + _ = primalityCache.TryAdd(composite, isPrime); } - } - - 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(); - - long timestamp = 0L; // reused for deadline checks. - 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; - } - - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - 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 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); + // else + // { + // Console.WriteLine($"isPrimeHits = {++_isPrimeHits}"); + // } + if (isPrime) { AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - } - else if (composite > 1UL) - { - cofactor = checked(cofactor * composite); - } - } - - ArrayPool pool = ThreadStaticPools.FactorEntryPool; - if (useDictionary) - { - if ((counts is null || counts.Count == 0) && cofactor == value) - { - return PartialFactorResult.Rent(null, value, false, 0); - } - - if (counts is null || counts.Count == 0) - { - return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - } - - FactorEntry[] factors = pool.Rent(counts.Count); - int index = 0; - foreach (KeyValuePair entry in counts) - { - factors[index] = new FactorEntry(entry.Key, entry.Value); - index++; - } - - Span span = factors.AsSpan(0, index); - span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - - bool fullyFactored = cofactor == 1UL; - return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); - } - - int actualCount = 0; - for (int i = 0; i < factorCount; i++) - { - if (primeSlots[i] != 0UL && exponentSlots[i] != 0) - { - actualCount++; - } - } - - if (actualCount == 0) - { - if (cofactor == value) - { - return PartialFactorResult.Rent(null, value, false, 0); - } - - return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - } - - FactorEntry[] array = pool.Rent(actualCount); - int arrayIndex = 0; - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeSlots[i]; - int exponentValue = exponentSlots[i]; - if (primeValue == 0UL || exponentValue == 0) - { continue; } - array[arrayIndex] = new FactorEntry(primeValue, exponentValue); - arrayIndex++; + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + 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 cofactor = 1UL; + int pendingCount = pending.Count; + for (int i = 0; i < pendingCount; i++) + { + ulong composite = pending[i]; + 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); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * composite); + } + } + + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + if (useDictionary) + { + if ((counts is null || counts.Count == 0) && cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; } - 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); + 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++; + } + + Span span = factors.AsSpan(0, index); + span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + + bool fullyFactored = cofactor == 1UL; + result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); + goto ReturnResult; } - finally + + int actualCount = 0; + for (int i = 0; i < factorCount; i++) { - pending.Clear(); - counts?.Clear(); - compositeStack?.Clear(); + if (primeSlots[i] != 0UL && exponentSlots[i] != 0) + { + actualCount++; + } } + + if (actualCount == 0) + { + if (cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; + } + + 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; + } + + 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; + result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); + + ReturnResult: + // ThreadStaticPools.UlongPool.Return(primeSlotsArray); + // ThreadStaticPools.IntPool.Return(exponentSlotsArray); + + if (counts is not null) + { + ThreadStaticPools.ReturnUlongIntDictionary(counts); + } + + if (compositeStack is not null) + { + compositeStack.Clear(); + ThreadStaticPools.ReturnUlongStack(compositeStack); + } + + pending.Clear(); + ThreadStaticPools.ReturnUlongList(pending); + 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( - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining) + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) { factorCount = 0; remaining = value; - if (value <= 1UL) - { - return true; - } + // value will never <= 1 in production code + // if (value <= 1UL) + // { + // return true; + // } int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); - if (capacity == 0) - { - return false; - } + // capacity will never equal 0 in production code + // if (capacity == 0) + // { + // return false; + // } uint[] primes = PrimesGenerator.SmallPrimes; ulong[] squares = PrimesGenerator.SmallPrimesPow2; int primeCount = primes.Length; ulong remainingLocal = value; uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + 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++) { @@ -1043,31 +1118,69 @@ internal static partial class PrimeOrderCalculator } ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remainingLocal) + // primeSquare will never == 0 in production code + // if (primeSquare != 0 && primeSquare > remainingLocal) + if (primeSquare > remainingLocal) { break; } - ulong primeValue = primeCandidate; - if ((remainingLocal % primeValue) != 0UL) + // primeCandidate will never equal 0 in production code + // if (primeCandidate == 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; } + // 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 { - remainingLocal /= primeValue; + remainingLocal /= primeCandidate; exponent++; } - while ((remainingLocal % primeValue) == 0UL); + while ((remainingLocal % primeCandidate) == 0UL); if (factorCount >= capacity) { - return false; + throw new InvalidOperationException($"Capacity is smaller than factor count"); } - primeTargets[factorCount] = primeValue; + primeTargets[factorCount] = primeCandidate; exponentTargets[factorCount] = exponent; factorCount++; } @@ -1118,14 +1231,18 @@ internal static partial class PrimeOrderCalculator private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) { - Span primeBuffer = stackalloc ulong[GpuSmallPrimeFactorSlots]; - Span exponentBuffer = stackalloc int[GpuSmallPrimeFactorSlots]; + 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; } @@ -1142,79 +1259,81 @@ internal static partial class PrimeOrderCalculator 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; - } + 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) @@ -1318,7 +1437,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..ed00820 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -79,6 +79,142 @@ 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; + + [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 fa0d030..c09ce8e 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; @@ -346,55 +345,55 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount); + + ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + Span 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; - } - - int windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } - - while (((exponent >> windowStart) & 1UL) == 0UL) - { - windowStart++; - } - - 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); - - index = windowStart - 1; - } - - if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold) - { - ArrayPool.Shared.Return(oddPowers); - } - - if (keepMontgomery) + while (index >= 0) { - return result; + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + 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); + + index = windowStart - 1; } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + ThreadStaticPools.UlongPool.Return(oddPowersArray); + + if (!keepMontgomery) + { + result = result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + return result; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -533,26 +532,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)]