From fa899587365415529cd591589bc45b5c92b53bb3 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Mon, 10 Nov 2025 14:30:30 +0100 Subject: [PATCH] Stabilization improvements. Stream will be released when releasing GPU lease to prevent application crashes. --- PerfectNumbers.Core/Gpu/AcceleratorPool.cs | 25 +- .../Accelerators/Pow2MontgomeryAccelerator.cs | 34 +- .../Gpu/Kernels/PrimeOrderKernels.cs | 2569 +++++++++-------- .../MersenneNumberDivisorByDivisorTester.cs | 2 +- PerfectNumbers.Core/PerfectNumberConstants.cs | 3 +- .../PrimeOrderCalculator.Cpu.cs | 20 +- .../PrimeOrderCalculator.Gpu.cs | 45 +- PerfectNumbers.Core/PrimeTester.cs | 10 +- 8 files changed, 1414 insertions(+), 1294 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/AcceleratorPool.cs b/PerfectNumbers.Core/Gpu/AcceleratorPool.cs index e6b1b5d..a796464 100644 --- a/PerfectNumbers.Core/Gpu/AcceleratorPool.cs +++ b/PerfectNumbers.Core/Gpu/AcceleratorPool.cs @@ -3,21 +3,21 @@ using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu { - public struct AcceleratorPool(int capacity) + public struct AcceleratorPool { - internal readonly Accelerator[] Accelerators = [.. Enumerable.Range(1, capacity).Select(_ => { - return SharedGpuContext.CreateAccelerator(); - })]; + internal readonly Context Context; + internal readonly Device Device; + internal readonly Accelerator[] Accelerators; private int _index; public Accelerator Rent() { - Atomic.CompareExchange(ref _index, capacity, 0); + Atomic.CompareExchange(ref _index, Capacity, 0); var index = Atomic.Add(ref _index, 1); - if (index >= capacity) + if (index >= Capacity) { - index -= capacity; + index -= Capacity; } return Accelerators[index]; @@ -29,5 +29,16 @@ namespace PerfectNumbers.Core.Gpu } public static readonly AcceleratorPool Shared = new(PerfectNumberConstants.RollingAccelerators); + private readonly int Capacity; + + public AcceleratorPool(int capacity) + { + Capacity = capacity; + var context = Context.CreateDefault(); + var device = context.GetPreferredDevice(false); + Context = context; + Device = device; + Accelerators = [.. Enumerable.Range(0, capacity).Select(_ => device.CreateAccelerator(context))]; + } } } \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs index eb3f6bf..1bd5cdb 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/Pow2MontgomeryAccelerator.cs @@ -8,7 +8,7 @@ namespace PerfectNumbers.Core.Gpu.Accelerators; public sealed class Pow2MontgomeryAccelerator { #region Pool - + [ThreadStatic] private static Queue? _pool; @@ -21,7 +21,11 @@ public sealed class Pow2MontgomeryAccelerator if (!pool.TryDequeue(out var gpu)) { Accelerator accelerator = _accelerators.Rent(); - gpu = new(accelerator); + gpu = new(accelerator, PerfectNumberConstants.DefaultFactorsBuffer); + } + else + { + gpu.Stream = gpu.Accelerator.CreateStream(); } return gpu; @@ -40,17 +44,33 @@ public sealed class Pow2MontgomeryAccelerator public readonly Accelerator Accelerator; public readonly MemoryBuffer1D Input; public readonly MemoryBuffer1D Output; + public MemoryBuffer1D, Stride1D.Dense> ToTestOnDevice; + public Dictionary ToTestOnHost = []; public AcceleratorStream? Stream; public readonly Action, MontgomeryDivisorData, ArrayView1D> ConvertToStandardKernel; public readonly Action, MontgomeryDivisorData, ArrayView1D> KeepMontgomeryKernel; - public Pow2MontgomeryAccelerator(Accelerator accelerator) + public readonly Kernel CheckFactorsKernel; + + public void EnsureCapacity(int factorsCount) + { + if (factorsCount > ToTestOnDevice.Length) + { + Accelerator accelerator = Accelerator; + ToTestOnDevice.Dispose(); + ToTestOnDevice = accelerator.Allocate1D>(factorsCount); + } + } + + public Pow2MontgomeryAccelerator(Accelerator accelerator, int factorsCount) { Accelerator = accelerator; + Stream = accelerator.CreateStream(); + Input = accelerator.Allocate1D(1); Output = accelerator.Allocate1D(1); - Stream = accelerator.CreateStream(); + ToTestOnDevice = accelerator.Allocate1D>(factorsCount); var keepKernel = accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery); var keepLauncher = KernelUtil.GetKernel(keepKernel); @@ -61,6 +81,9 @@ public sealed class Pow2MontgomeryAccelerator var convertLauncher = KernelUtil.GetKernel(convertKernel); ConvertToStandardKernel = convertLauncher.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); + + var checkFactorsKernel = accelerator.LoadStreamKernel, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D>(PrimeOrderGpuHeuristics.CheckFactorsKernel); + CheckFactorsKernel = KernelUtil.GetKernel(checkFactorsKernel); } public void Dispose() @@ -68,7 +91,8 @@ public sealed class Pow2MontgomeryAccelerator Input.Dispose(); Output.Dispose(); Stream?.Dispose(); + ToTestOnDevice.Dispose(); // Accelerator is shared with other threads // Accelerator.Dispose(); } -} +} \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs index 5c5ae14..3ba4077 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs @@ -7,144 +7,183 @@ namespace PerfectNumbers.Core.Gpu; internal static partial class PrimeOrderGpuHeuristics { + internal static void CheckFactorsKernel( + int count, + ulong phi, + ArrayView1D, Stride1D.Dense> toTest, + MontgomeryDivisorData divisorData, + ArrayView1D result + ) + { + ulong order = phi; + int entryCount = count; + + for (int i = 0; i < entryCount; i++) + { + (ulong primeFactor, int exponent) = toTest[i]; + + for (int iteration = 0; iteration < exponent; iteration++) + { + ulong candidate = order / primeFactor; + ulong remainder = order - (candidate * primeFactor); + // ulong remainder = order % primeFactor; + if (remainder != 0UL) + { + break; + } + + // ulong candidate = order / primeFactor; + if (ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisorData, candidate) == 1UL) + { + order = candidate; + continue; + } + + break; + } + } + + result[0] = order; + } + /// This kernel always sets the result and operates within the allowed / required bounds, never accessing elements outside of them. /// It sets the value only when it needs to. The callers don't need to clear the output / input buffers. - internal static void PartialFactorKernel( - Index1D index, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - int slotCount, - ulong value, - uint limit, - ArrayView1D factorsOut, - ArrayView1D exponentsOut, - ArrayView1D countOut, - ArrayView1D remainingOut, - ArrayView1D fullyFactoredOut) - { - if (index != 0) - { - return; - } + internal static void PartialFactorKernel( + Index1D index, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + int slotCount, + ulong value, + uint limit, + ArrayView1D factorsOut, + ArrayView1D exponentsOut, + ArrayView1D countOut, + ArrayView1D remainingOut, + ArrayView1D fullyFactoredOut) + { + if (index != 0) + { + return; + } - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - ulong remainingLocal = value; - int count = 0; + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + ulong remainingLocal = value; + int count = 0; - for (int i = 0; i < primeCount && count < slotCount; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > effectiveLimit) - { - break; - } + for (int i = 0; i < primeCount && count < slotCount; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remainingLocal) - { - break; - } + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remainingLocal) + { + break; + } - ulong primeValue = primeCandidate; - if (primeValue == 0UL || (remainingLocal % primeValue) != 0UL) - { - continue; - } + ulong primeValue = primeCandidate; + if (primeValue == 0UL || (remainingLocal % primeValue) != 0UL) + { + continue; + } - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); - factorsOut[count] = primeValue; - exponentsOut[count] = exponent; - count++; - } + factorsOut[count] = primeValue; + exponentsOut[count] = exponent; + count++; + } - countOut[0] = count; - remainingOut[0] = remainingLocal; - fullyFactoredOut[0] = remainingLocal == 1UL ? (byte)1 : (byte)0; - } + countOut[0] = count; + remainingOut[0] = remainingLocal; + fullyFactoredOut[0] = remainingLocal == 1UL ? (byte)1 : (byte)0; + } - internal enum PrimeOrderKernelStatus : byte - { - Fallback = 0, - Found = 1, - HeuristicUnresolved = 2, - PollardOverflow = 3, - FactoringFailure = 4, - } + internal enum PrimeOrderKernelStatus : byte + { + Fallback = 0, + Found = 1, + HeuristicUnresolved = 2, + PollardOverflow = 3, + FactoringFailure = 4, + } - private readonly struct CandidateKey - { - public CandidateKey(int primary, long secondary, long tertiary) - { - Primary = primary; - Secondary = secondary; - Tertiary = tertiary; - } + private readonly struct CandidateKey + { + public CandidateKey(int primary, long secondary, long tertiary) + { + Primary = primary; + Secondary = secondary; + Tertiary = tertiary; + } - public int Primary { get; } + public int Primary { get; } - public long Secondary { get; } + public long Secondary { get; } - public long Tertiary { get; } + public long Tertiary { get; } - public int CompareTo(CandidateKey other) - { - if (Primary != other.Primary) - { - return Primary.CompareTo(other.Primary); - } + public int CompareTo(CandidateKey other) + { + if (Primary != other.Primary) + { + return Primary.CompareTo(other.Primary); + } - if (Secondary != other.Secondary) - { - return Secondary.CompareTo(other.Secondary); - } + if (Secondary != other.Secondary) + { + return Secondary.CompareTo(other.Secondary); + } - return Tertiary.CompareTo(other.Tertiary); - } - } + return Tertiary.CompareTo(other.Tertiary); + } + } /// This kernel doesn't always set the result of the first element. It sets the value only when it needs to. The callers /// must clean the output buffer before the call to get a deterministic result. - internal static void CalculateOrderKernel( - Index1D index, - ulong prime, - OrderKernelConfig config, - MontgomeryDivisorData divisor, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - OrderKernelBuffers buffers) - { - if (index != 0) - { - return; - } + internal static void CalculateOrderKernel( + Index1D index, + ulong prime, + OrderKernelConfig config, + MontgomeryDivisorData divisor, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + OrderKernelBuffers buffers) + { + if (index != 0) + { + return; + } // TODO: Modify this kernel to always set the output result so that callers never need to clear any buffers and remove // any existing cleaning of the buffers after that. ArrayView1D phiFactors = buffers.PhiFactors; - ArrayView1D phiExponents = buffers.PhiExponents; - ArrayView1D workFactors = buffers.WorkFactors; - ArrayView1D workExponents = buffers.WorkExponents; - ArrayView1D candidates = buffers.Candidates; - ArrayView1D stackIndex = buffers.StackIndex; - ArrayView1D stackExponent = buffers.StackExponent; - ArrayView1D stackProduct = buffers.StackProduct; - ArrayView1D resultOut = buffers.Result; - ArrayView1D statusOut = buffers.Status; + ArrayView1D phiExponents = buffers.PhiExponents; + ArrayView1D workFactors = buffers.WorkFactors; + ArrayView1D workExponents = buffers.WorkExponents; + ArrayView1D candidates = buffers.Candidates; + ArrayView1D stackIndex = buffers.StackIndex; + ArrayView1D stackExponent = buffers.StackExponent; + ArrayView1D stackProduct = buffers.StackProduct; + ArrayView1D resultOut = buffers.Result; + ArrayView1D statusOut = buffers.Status; - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong previousOrder = config.PreviousOrder; - byte hasPreviousOrder = config.HasPreviousOrder; - int maxPowChecks = config.MaxPowChecks; - int mode = config.Mode; + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong previousOrder = config.PreviousOrder; + byte hasPreviousOrder = config.HasPreviousOrder; + int maxPowChecks = config.MaxPowChecks; + int mode = config.Mode; statusOut[0] = (byte)PrimeOrderKernelStatus.Fallback; @@ -159,1005 +198,1005 @@ internal static partial class PrimeOrderGpuHeuristics return; } - ulong phi = prime - 1UL; - - int phiFactorCount = FactorWithSmallPrimes(phi, limit, primes, squares, primeCount, phiFactors, phiExponents, out ulong phiRemaining); - if (phiRemaining != 1UL) - { - // Reuse stackProduct as the Pollard Rho stack so factoring stays within this kernel. - if (!TryFactorWithPollardKernel( - phiRemaining, - limit, - primes, - squares, - primeCount, - phiFactors, - phiExponents, - ref phiFactorCount, - stackProduct, - statusOut)) - { - if (statusOut[0] != (byte)PrimeOrderKernelStatus.PollardOverflow) - { - resultOut[0] = CalculateByDoublingKernel(prime); - } - - return; - } - } - - SortFactors(phiFactors, phiExponents, phiFactorCount); - - if (TrySpecialMaxKernel(phi, prime, phiFactors, phiFactorCount, divisor)) - { - resultOut[0] = phi; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - ulong candidateOrder = InitializeStartingOrderKernel(prime, phi, divisor); - candidateOrder = ExponentLoweringKernel(candidateOrder, prime, phiFactors, phiExponents, phiFactorCount, divisor); - - if (TryConfirmOrderKernel( - prime, - candidateOrder, - divisor, - limit, - primes, - squares, - primeCount, - workFactors, - workExponents, - stackProduct, - statusOut)) - { - resultOut[0] = candidateOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - bool isStrict = mode == 1; - if (isStrict) - { - ulong strictOrder = CalculateByDoublingKernel(prime); - resultOut[0] = strictOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - if (TryHeuristicFinishKernel( - prime, - candidateOrder, - previousOrder, - hasPreviousOrder, - divisor, - limit, - maxPowChecks, - primes, - squares, - primeCount, - workFactors, - workExponents, - candidates, - stackIndex, - stackExponent, - stackProduct, - statusOut, - out ulong confirmedOrder)) - { - resultOut[0] = confirmedOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - ulong fallbackOrder = CalculateByDoublingKernel(prime); - resultOut[0] = fallbackOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.HeuristicUnresolved; - } - - private static int FactorWithSmallPrimes( - ulong value, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - out ulong remaining) - { - remaining = value; - int factorCount = 0; - long factorLength = factors.Length; - long exponentLength = exponents.Length; - int capacity = factorLength < exponentLength ? (int)factorLength : (int)exponentLength; - - for (int i = 0; i < primeCount && remaining > 1UL && factorCount < capacity; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate == 0U || primeCandidate > limit) - { - break; - } - - ulong square = squares[i]; - if (square != 0UL && square > remaining) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == 0UL); - - factors[factorCount] = primeValue; - exponents[factorCount] = exponent; - factorCount++; - } - - for (int i = factorCount; i < factors.Length; i++) - { - factors[i] = 0UL; - } - - for (int i = factorCount; i < exponents.Length; i++) - { - exponents[i] = 0; - } - - return factorCount; - } - - private static void SortFactors(ArrayView1D factors, ArrayView1D exponents, int count) - { - for (int i = 1; i < count; i++) - { - ulong factor = factors[i]; - int exponent = exponents[i]; - int j = i - 1; - - while (j >= 0 && factors[j] > factor) - { - factors[j + 1] = factors[j]; - exponents[j + 1] = exponents[j]; - j--; - } - - factors[j + 1] = factor; - exponents[j + 1] = exponent; - } - } - - private static bool TrySpecialMaxKernel( - ulong phi, - ulong prime, - ArrayView1D factors, - int factorCount, - in MontgomeryDivisorData divisor) - { - for (int i = 0; i < factorCount; i++) - { - ulong factor = factors[i]; - if (factor <= 1UL) - { - continue; - } - - ulong reduced = phi / factor; - if (ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, reduced) == 1UL) - { - return false; - } - } - - return true; - } - - private static ulong InitializeStartingOrderKernel(ulong prime, ulong phi, in MontgomeryDivisorData divisor) - { - ulong order = phi; - ulong residue = prime & 7UL; - if (residue == 1UL || residue == 7UL) - { - ulong half = phi >> 1; - if (Pow2EqualsOneKernel(half, divisor)) - { - order = half; - } - } - - return order; - } - - private static ulong ExponentLoweringKernel( - ulong order, - ulong prime, - ArrayView1D factors, - ArrayView1D exponents, - int factorCount, - in MontgomeryDivisorData divisor) - { - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factors[i]; - int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } - - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } - - ulong reduced = order / primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) - { - order = reduced; - continue; - } - - break; - } - } - - return order; - } - - private static bool TryConfirmOrderKernel( - ulong prime, - ulong order, - in MontgomeryDivisorData divisor, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ArrayView1D compositeStack, - ArrayView1D statusOut) - { - if (order == 0UL) - { - return false; - } - - if (!Pow2EqualsOneKernel(order, divisor)) - { - return false; - } - - int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); - if (remaining != 1UL) - { - if (!TryFactorWithPollardKernel( - remaining, - limit, - primes, - squares, - primeCount, - factors, - exponents, - ref factorCount, - compositeStack, - statusOut)) - { - return false; - } - } - - SortFactors(factors, exponents, factorCount); - - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factors[i]; - int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } - - ulong reduced = order; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) - { - return false; - } - } - } - - return true; - } - - private static bool TryHeuristicFinishKernel( - ulong prime, - ulong order, - ulong previousOrder, - byte hasPreviousOrder, - in MontgomeryDivisorData divisor, - uint limit, - int maxPowChecks, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D workFactors, - ArrayView1D workExponents, - ArrayView1D candidates, - ArrayView1D stackIndex, - ArrayView1D stackExponent, - ArrayView1D stackProduct, - ArrayView1D statusOut, - out ulong confirmedOrder) - { - confirmedOrder = 0UL; - - if (order <= 1UL) - { - return false; - } - - int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, workFactors, workExponents, out ulong remaining); - if (remaining != 1UL) - { - // Reuse stackProduct as the Pollard Rho stack while factoring the order candidates. - if (!TryFactorWithPollardKernel( - remaining, - limit, - primes, - squares, - primeCount, - workFactors, - workExponents, - ref factorCount, - stackProduct, - statusOut)) - { - return false; - } - } - - SortFactors(workFactors, workExponents, factorCount); - - long candidateCapacity = candidates.Length; - int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit; - int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit); - if (candidateCount == 0) - { - return false; - } - - SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); - - int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks; - if (powBudget <= 0) - { - powBudget = candidateCount; - } - - int powUsed = 0; - - for (int i = 0; i < candidateCount && powUsed < powBudget; i++) - { - ulong candidate = candidates[i]; - if (candidate <= 1UL) - { - continue; - } - - if (powUsed >= powBudget) - { - break; - } - - powUsed++; - if (!Pow2EqualsOneKernel(candidate, divisor)) - { - continue; - } - - if (!TryConfirmCandidateKernel( - prime, - candidate, - divisor, - limit, - primes, - squares, - primeCount, - workFactors, - workExponents, - stackProduct, - statusOut, - ref powUsed, - powBudget)) - { - continue; - } - - confirmedOrder = candidate; - return true; - } - - return false; - } - - private static int BuildCandidatesKernel( - ulong order, - ArrayView1D factors, - ArrayView1D exponents, - int factorCount, - ArrayView1D candidates, - ArrayView1D stackIndex, - ArrayView1D stackExponent, - ArrayView1D stackProduct, - int limit) - { - long stackIndexLength = stackIndex.Length; - long stackExponentLength = stackExponent.Length; - long stackProductLength = stackProduct.Length; - if (factorCount == 0 || limit <= 0 || stackIndexLength == 0L || stackExponentLength == 0L || stackProductLength == 0L) - { - return 0; - } - - int stackCapacity = stackIndexLength < stackExponentLength ? (int)stackIndexLength : (int)stackExponentLength; - if (stackProductLength < stackCapacity) - { - stackCapacity = (int)stackProductLength; - } - - int candidateCount = 0; - int stackTop = 0; - - stackIndex[0] = 0; - stackExponent[0] = 0; - stackProduct[0] = 1UL; - stackTop = 1; - - while (stackTop > 0) - { - stackTop--; - int index = stackIndex[stackTop]; - int exponent = stackExponent[stackTop]; - ulong product = stackProduct[stackTop]; - - if (index >= factorCount) - { - if (product != 1UL && product != order && candidateCount < limit) - { - ulong candidate = order / product; - if (candidate > 1UL && candidate < order) - { - candidates[candidateCount] = candidate; - candidateCount++; - } - } - - continue; - } - - int maxExponent = exponents[index]; - if (exponent > maxExponent) - { - continue; - } - - if (stackTop >= stackCapacity) - { - return candidateCount; - } - - stackIndex[stackTop] = index + 1; - stackExponent[stackTop] = 0; - stackProduct[stackTop] = product; - stackTop++; - - if (exponent == maxExponent) - { - continue; - } - - ulong primeFactor = factors[index]; - if (primeFactor == 0UL || product > order / primeFactor) - { - continue; - } - - if (stackTop >= stackCapacity) - { - return candidateCount; - } - - stackIndex[stackTop] = index; - stackExponent[stackTop] = exponent + 1; - stackProduct[stackTop] = product * primeFactor; - stackTop++; - } - - return candidateCount; - } - - private static void SortCandidatesKernel( - ulong prime, - ulong previousOrder, - bool hasPrevious, - ArrayView1D candidates, - int count) - { - for (int i = 1; i < count; i++) - { - ulong value = candidates[i]; - CandidateKey key = BuildCandidateKey(value, prime, previousOrder, hasPrevious); - int j = i - 1; - - while (j >= 0) - { - CandidateKey other = BuildCandidateKey(candidates[j], prime, previousOrder, hasPrevious); - if (other.CompareTo(key) <= 0) - { - break; - } - - candidates[j + 1] = candidates[j]; - j--; - } - - candidates[j + 1] = value; - } - } - - private static CandidateKey BuildCandidateKey(ulong value, ulong prime, ulong previousOrder, bool hasPrevious) - { - int group = GetGroup(value, prime); - if (group == 0) - { - return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); - } - - ulong reference = hasPrevious ? previousOrder : 0UL; - bool isGe = !hasPrevious || value >= reference; - int previousGroup = hasPrevious ? GetGroup(reference, prime) : 1; - int primary = ComputePrimary(group, isGe, previousGroup); - long secondary; - long tertiary; - - if (group == 3) - { - secondary = -(long)value; - tertiary = -(long)value; - } - else - { - ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; - secondary = (long)distance; - tertiary = (long)value; - } - - return new CandidateKey(primary, secondary, tertiary); - } - - private static int GetGroup(ulong value, ulong prime) - { - ulong threshold1 = prime >> 3; - if (value <= threshold1) - { - return 1; - } - - ulong threshold2 = prime >> 2; - if (value <= threshold2) - { - return 2; - } - - ulong threshold3 = (prime * 3UL) >> 3; - if (value <= threshold3) - { - return 3; - } - - return 0; - } - - private static int ComputePrimary(int group, bool isGe, int previousGroup) - { - int groupOffset; - switch (group) - { - case 1: - groupOffset = 0; - break; - case 2: - groupOffset = 2; - break; - case 3: - groupOffset = 4; - break; - default: - groupOffset = 6; - break; - } - - if (group == previousGroup) - { - if (group == 3) - { - return groupOffset + (isGe ? 0 : 3); - } - - return groupOffset + (isGe ? 0 : 1); - } - - return groupOffset + (isGe ? 0 : 1); - } - - private static bool TryConfirmCandidateKernel( - ulong prime, - ulong candidate, - in MontgomeryDivisorData divisor, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ArrayView1D compositeStack, - ArrayView1D statusOut, - ref int powUsed, - int powBudget) - { - int factorCount = FactorWithSmallPrimes(candidate, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); - if (remaining != 1UL) - { - if (!TryFactorWithPollardKernel( - remaining, - limit, - primes, - squares, - primeCount, - factors, - exponents, - ref factorCount, - compositeStack, - statusOut)) - { - return false; - } - } - - SortFactors(factors, exponents, factorCount); - - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factors[i]; - int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } - - ulong reduced = candidate; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (powBudget > 0 && powUsed >= powBudget) - { - return false; - } - - powUsed++; - if (Pow2EqualsOneKernel(reduced, divisor)) - { - return false; - } - } - } - - return true; - } - - private static bool TryFactorWithPollardKernel( - ulong initial, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ref int factorCount, - ArrayView1D compositeStack, - ArrayView1D statusOut) - { - if (initial <= 1UL) - { - return true; - } - - int stackCapacity = (int)compositeStack.Length; - if (stackCapacity <= 0) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - int stackTop = 0; - compositeStack[stackTop] = initial; - stackTop++; - - while (stackTop > 0) - { - stackTop--; - ulong composite = compositeStack[stackTop]; - if (composite <= 1UL) - { - continue; - } - - if (!PeelSmallPrimesKernel( - composite, - limit, - primes, - squares, - primeCount, - factors, - exponents, - ref factorCount, - out ulong remaining)) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - if (remaining == 1UL) - { - continue; - } - - ulong factor = PollardRhoKernel(remaining); - if (factor <= 1UL || factor == remaining) - { - if (!TryAppendFactorKernel(factors, exponents, ref factorCount, remaining, 1)) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - continue; - } - - ulong quotient = remaining / factor; - if (stackTop + 2 > stackCapacity) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - compositeStack[stackTop] = factor; - compositeStack[stackTop + 1] = quotient; - stackTop += 2; - } - - return true; - } - - private static bool PeelSmallPrimesKernel( - ulong value, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ref int factorCount, - out ulong remaining) - { - ulong remainingLocal = value; - - for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate == 0U || primeCandidate > limit) - { - break; - } - - ulong square = squares[i]; - if (square != 0UL && square > remainingLocal) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remainingLocal % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); - - if (!TryAppendFactorKernel(factors, exponents, ref factorCount, primeValue, exponent)) - { - remaining = value; - return false; - } - } - - remaining = remainingLocal; - return true; - } - - private static bool TryAppendFactorKernel( - ArrayView1D factors, - ArrayView1D exponents, - ref int count, - ulong prime, - int exponent) - { - if (prime <= 1UL || exponent <= 0) - { - return true; - } - - for (int i = 0; i < count; i++) - { - if (factors[i] == prime) - { - exponents[i] += exponent; - return true; - } - } - - int capacity = (int)factors.Length; - if (count >= capacity) - { - return false; - } - - factors[count] = prime; - exponents[count] = exponent; - count++; - return true; - } - - private static ulong MulModKernel(ulong left, ulong right, ulong modulus) - { - GpuUInt128 product = new GpuUInt128(left); - return product.MulMod(right, modulus); - } - - private static ulong PollardRhoKernel(ulong value) - { - if ((value & 1UL) == 0UL) - { - return 2UL; - } - - ulong c = 1UL; - while (true) - { - ulong x = 2UL; - ulong y = 2UL; - ulong d = 1UL; - - while (d == 1UL) - { - x = AdvancePolynomialKernel(x, c, value); - y = AdvancePolynomialKernel(y, c, value); - y = AdvancePolynomialKernel(y, c, value); - - ulong diff = x > y ? x - y : y - x; - d = BinaryGcdKernel(diff, value); - } - - if (d == value) - { - c++; - if (c == 0UL) - { - c = 1UL; - } - - continue; - } - - return d; - } - } - - private static ulong AdvancePolynomialKernel(ulong x, ulong c, ulong modulus) - { - ulong squared = MulModKernel(x, x, modulus); - GpuUInt128 accumulator = new GpuUInt128(squared); - accumulator.AddMod(c, modulus); - return accumulator.Low; - } - - private static ulong BinaryGcdKernel(ulong a, ulong b) - { - if (a == 0UL) - { - return b; - } - - if (b == 0UL) - { - return a; - } - - int shift = BitOperations.TrailingZeroCount(a | b); - ulong aLocal = a >> BitOperations.TrailingZeroCount(a); - ulong bLocal = b; - - while (true) - { - bLocal >>= BitOperations.TrailingZeroCount(bLocal); - if (aLocal > bLocal) - { - ulong temp = aLocal; - aLocal = bLocal; - bLocal = temp; - } - - bLocal -= aLocal; - if (bLocal == 0UL) - { - return aLocal << shift; - } - } - } - - private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) - { - return ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent) == 1UL; - } - - private static ulong CalculateByDoublingKernel(ulong prime) - { - ulong order = 1UL; - ulong value = 2UL % prime; - - while (value != 1UL) - { - value <<= 1; - if (value >= prime) - { - value -= prime; - } - - order++; - } - - return order; - } + ulong phi = prime - 1UL; + + int phiFactorCount = FactorWithSmallPrimes(phi, limit, primes, squares, primeCount, phiFactors, phiExponents, out ulong phiRemaining); + if (phiRemaining != 1UL) + { + // Reuse stackProduct as the Pollard Rho stack so factoring stays within this kernel. + if (!TryFactorWithPollardKernel( + phiRemaining, + limit, + primes, + squares, + primeCount, + phiFactors, + phiExponents, + ref phiFactorCount, + stackProduct, + statusOut)) + { + if (statusOut[0] != (byte)PrimeOrderKernelStatus.PollardOverflow) + { + resultOut[0] = CalculateByDoublingKernel(prime); + } + + return; + } + } + + SortFactors(phiFactors, phiExponents, phiFactorCount); + + if (TrySpecialMaxKernel(phi, prime, phiFactors, phiFactorCount, divisor)) + { + resultOut[0] = phi; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + ulong candidateOrder = InitializeStartingOrderKernel(prime, phi, divisor); + candidateOrder = ExponentLoweringKernel(candidateOrder, prime, phiFactors, phiExponents, phiFactorCount, divisor); + + if (TryConfirmOrderKernel( + prime, + candidateOrder, + divisor, + limit, + primes, + squares, + primeCount, + workFactors, + workExponents, + stackProduct, + statusOut)) + { + resultOut[0] = candidateOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + bool isStrict = mode == 1; + if (isStrict) + { + ulong strictOrder = CalculateByDoublingKernel(prime); + resultOut[0] = strictOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + if (TryHeuristicFinishKernel( + prime, + candidateOrder, + previousOrder, + hasPreviousOrder, + divisor, + limit, + maxPowChecks, + primes, + squares, + primeCount, + workFactors, + workExponents, + candidates, + stackIndex, + stackExponent, + stackProduct, + statusOut, + out ulong confirmedOrder)) + { + resultOut[0] = confirmedOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + ulong fallbackOrder = CalculateByDoublingKernel(prime); + resultOut[0] = fallbackOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.HeuristicUnresolved; + } + + private static int FactorWithSmallPrimes( + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + out ulong remaining) + { + remaining = value; + int factorCount = 0; + long factorLength = factors.Length; + long exponentLength = exponents.Length; + int capacity = factorLength < exponentLength ? (int)factorLength : (int)exponentLength; + + for (int i = 0; i < primeCount && remaining > 1UL && factorCount < capacity; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > limit) + { + break; + } + + ulong square = squares[i]; + if (square != 0UL && square > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + factors[factorCount] = primeValue; + exponents[factorCount] = exponent; + factorCount++; + } + + for (int i = factorCount; i < factors.Length; i++) + { + factors[i] = 0UL; + } + + for (int i = factorCount; i < exponents.Length; i++) + { + exponents[i] = 0; + } + + return factorCount; + } + + private static void SortFactors(ArrayView1D factors, ArrayView1D exponents, int count) + { + for (int i = 1; i < count; i++) + { + ulong factor = factors[i]; + int exponent = exponents[i]; + int j = i - 1; + + while (j >= 0 && factors[j] > factor) + { + factors[j + 1] = factors[j]; + exponents[j + 1] = exponents[j]; + j--; + } + + factors[j + 1] = factor; + exponents[j + 1] = exponent; + } + } + + private static bool TrySpecialMaxKernel( + ulong phi, + ulong prime, + ArrayView1D factors, + int factorCount, + in MontgomeryDivisorData divisor) + { + for (int i = 0; i < factorCount; i++) + { + ulong factor = factors[i]; + if (factor <= 1UL) + { + continue; + } + + ulong reduced = phi / factor; + if (ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, reduced) == 1UL) + { + return false; + } + } + + return true; + } + + private static ulong InitializeStartingOrderKernel(ulong prime, ulong phi, in MontgomeryDivisorData divisor) + { + ulong order = phi; + ulong residue = prime & 7UL; + if (residue == 1UL || residue == 7UL) + { + ulong half = phi >> 1; + if (Pow2EqualsOneKernel(half, divisor)) + { + order = half; + } + } + + return order; + } + + private static ulong ExponentLoweringKernel( + ulong order, + ulong prime, + ArrayView1D factors, + ArrayView1D exponents, + int factorCount, + in MontgomeryDivisorData divisor) + { + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factors[i]; + int exponent = exponents[i]; + if (primeFactor <= 1UL) + { + continue; + } + + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong reduced = order / primeFactor; + if (Pow2EqualsOneKernel(reduced, divisor)) + { + order = reduced; + continue; + } + + break; + } + } + + return order; + } + + private static bool TryConfirmOrderKernel( + ulong prime, + ulong order, + in MontgomeryDivisorData divisor, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ArrayView1D compositeStack, + ArrayView1D statusOut) + { + if (order == 0UL) + { + return false; + } + + if (!Pow2EqualsOneKernel(order, divisor)) + { + return false; + } + + int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); + if (remaining != 1UL) + { + if (!TryFactorWithPollardKernel( + remaining, + limit, + primes, + squares, + primeCount, + factors, + exponents, + ref factorCount, + compositeStack, + statusOut)) + { + return false; + } + } + + SortFactors(factors, exponents, factorCount); + + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factors[i]; + int exponent = exponents[i]; + if (primeFactor <= 1UL) + { + continue; + } + + ulong reduced = order; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (Pow2EqualsOneKernel(reduced, divisor)) + { + return false; + } + } + } + + return true; + } + + private static bool TryHeuristicFinishKernel( + ulong prime, + ulong order, + ulong previousOrder, + byte hasPreviousOrder, + in MontgomeryDivisorData divisor, + uint limit, + int maxPowChecks, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D workFactors, + ArrayView1D workExponents, + ArrayView1D candidates, + ArrayView1D stackIndex, + ArrayView1D stackExponent, + ArrayView1D stackProduct, + ArrayView1D statusOut, + out ulong confirmedOrder) + { + confirmedOrder = 0UL; + + if (order <= 1UL) + { + return false; + } + + int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, workFactors, workExponents, out ulong remaining); + if (remaining != 1UL) + { + // Reuse stackProduct as the Pollard Rho stack while factoring the order candidates. + if (!TryFactorWithPollardKernel( + remaining, + limit, + primes, + squares, + primeCount, + workFactors, + workExponents, + ref factorCount, + stackProduct, + statusOut)) + { + return false; + } + } + + SortFactors(workFactors, workExponents, factorCount); + + long candidateCapacity = candidates.Length; + int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit; + int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit); + if (candidateCount == 0) + { + return false; + } + + SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); + + int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks; + if (powBudget <= 0) + { + powBudget = candidateCount; + } + + int powUsed = 0; + + for (int i = 0; i < candidateCount && powUsed < powBudget; i++) + { + ulong candidate = candidates[i]; + if (candidate <= 1UL) + { + continue; + } + + if (powUsed >= powBudget) + { + break; + } + + powUsed++; + if (!Pow2EqualsOneKernel(candidate, divisor)) + { + continue; + } + + if (!TryConfirmCandidateKernel( + prime, + candidate, + divisor, + limit, + primes, + squares, + primeCount, + workFactors, + workExponents, + stackProduct, + statusOut, + ref powUsed, + powBudget)) + { + continue; + } + + confirmedOrder = candidate; + return true; + } + + return false; + } + + private static int BuildCandidatesKernel( + ulong order, + ArrayView1D factors, + ArrayView1D exponents, + int factorCount, + ArrayView1D candidates, + ArrayView1D stackIndex, + ArrayView1D stackExponent, + ArrayView1D stackProduct, + int limit) + { + long stackIndexLength = stackIndex.Length; + long stackExponentLength = stackExponent.Length; + long stackProductLength = stackProduct.Length; + if (factorCount == 0 || limit <= 0 || stackIndexLength == 0L || stackExponentLength == 0L || stackProductLength == 0L) + { + return 0; + } + + int stackCapacity = stackIndexLength < stackExponentLength ? (int)stackIndexLength : (int)stackExponentLength; + if (stackProductLength < stackCapacity) + { + stackCapacity = (int)stackProductLength; + } + + int candidateCount = 0; + int stackTop = 0; + + stackIndex[0] = 0; + stackExponent[0] = 0; + stackProduct[0] = 1UL; + stackTop = 1; + + while (stackTop > 0) + { + stackTop--; + int index = stackIndex[stackTop]; + int exponent = stackExponent[stackTop]; + ulong product = stackProduct[stackTop]; + + if (index >= factorCount) + { + if (product != 1UL && product != order && candidateCount < limit) + { + ulong candidate = order / product; + if (candidate > 1UL && candidate < order) + { + candidates[candidateCount] = candidate; + candidateCount++; + } + } + + continue; + } + + int maxExponent = exponents[index]; + if (exponent > maxExponent) + { + continue; + } + + if (stackTop >= stackCapacity) + { + return candidateCount; + } + + stackIndex[stackTop] = index + 1; + stackExponent[stackTop] = 0; + stackProduct[stackTop] = product; + stackTop++; + + if (exponent == maxExponent) + { + continue; + } + + ulong primeFactor = factors[index]; + if (primeFactor == 0UL || product > order / primeFactor) + { + continue; + } + + if (stackTop >= stackCapacity) + { + return candidateCount; + } + + stackIndex[stackTop] = index; + stackExponent[stackTop] = exponent + 1; + stackProduct[stackTop] = product * primeFactor; + stackTop++; + } + + return candidateCount; + } + + private static void SortCandidatesKernel( + ulong prime, + ulong previousOrder, + bool hasPrevious, + ArrayView1D candidates, + int count) + { + for (int i = 1; i < count; i++) + { + ulong value = candidates[i]; + CandidateKey key = BuildCandidateKey(value, prime, previousOrder, hasPrevious); + int j = i - 1; + + while (j >= 0) + { + CandidateKey other = BuildCandidateKey(candidates[j], prime, previousOrder, hasPrevious); + if (other.CompareTo(key) <= 0) + { + break; + } + + candidates[j + 1] = candidates[j]; + j--; + } + + candidates[j + 1] = value; + } + } + + private static CandidateKey BuildCandidateKey(ulong value, ulong prime, ulong previousOrder, bool hasPrevious) + { + int group = GetGroup(value, prime); + if (group == 0) + { + return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); + } + + ulong reference = hasPrevious ? previousOrder : 0UL; + bool isGe = !hasPrevious || value >= reference; + int previousGroup = hasPrevious ? GetGroup(reference, prime) : 1; + int primary = ComputePrimary(group, isGe, previousGroup); + long secondary; + long tertiary; + + if (group == 3) + { + secondary = -(long)value; + tertiary = -(long)value; + } + else + { + ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; + secondary = (long)distance; + tertiary = (long)value; + } + + return new CandidateKey(primary, secondary, tertiary); + } + + private static int GetGroup(ulong value, ulong prime) + { + ulong threshold1 = prime >> 3; + if (value <= threshold1) + { + return 1; + } + + ulong threshold2 = prime >> 2; + if (value <= threshold2) + { + return 2; + } + + ulong threshold3 = (prime * 3UL) >> 3; + if (value <= threshold3) + { + return 3; + } + + return 0; + } + + private static int ComputePrimary(int group, bool isGe, int previousGroup) + { + int groupOffset; + switch (group) + { + case 1: + groupOffset = 0; + break; + case 2: + groupOffset = 2; + break; + case 3: + groupOffset = 4; + break; + default: + groupOffset = 6; + break; + } + + if (group == previousGroup) + { + if (group == 3) + { + return groupOffset + (isGe ? 0 : 3); + } + + return groupOffset + (isGe ? 0 : 1); + } + + return groupOffset + (isGe ? 0 : 1); + } + + private static bool TryConfirmCandidateKernel( + ulong prime, + ulong candidate, + in MontgomeryDivisorData divisor, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ArrayView1D compositeStack, + ArrayView1D statusOut, + ref int powUsed, + int powBudget) + { + int factorCount = FactorWithSmallPrimes(candidate, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); + if (remaining != 1UL) + { + if (!TryFactorWithPollardKernel( + remaining, + limit, + primes, + squares, + primeCount, + factors, + exponents, + ref factorCount, + compositeStack, + statusOut)) + { + return false; + } + } + + SortFactors(factors, exponents, factorCount); + + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factors[i]; + int exponent = exponents[i]; + if (primeFactor <= 1UL) + { + continue; + } + + ulong reduced = candidate; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (powBudget > 0 && powUsed >= powBudget) + { + return false; + } + + powUsed++; + if (Pow2EqualsOneKernel(reduced, divisor)) + { + return false; + } + } + } + + return true; + } + + private static bool TryFactorWithPollardKernel( + ulong initial, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ref int factorCount, + ArrayView1D compositeStack, + ArrayView1D statusOut) + { + if (initial <= 1UL) + { + return true; + } + + int stackCapacity = (int)compositeStack.Length; + if (stackCapacity <= 0) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + int stackTop = 0; + compositeStack[stackTop] = initial; + stackTop++; + + while (stackTop > 0) + { + stackTop--; + ulong composite = compositeStack[stackTop]; + if (composite <= 1UL) + { + continue; + } + + if (!PeelSmallPrimesKernel( + composite, + limit, + primes, + squares, + primeCount, + factors, + exponents, + ref factorCount, + out ulong remaining)) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + if (remaining == 1UL) + { + continue; + } + + ulong factor = PollardRhoKernel(remaining); + if (factor <= 1UL || factor == remaining) + { + if (!TryAppendFactorKernel(factors, exponents, ref factorCount, remaining, 1)) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + continue; + } + + ulong quotient = remaining / factor; + if (stackTop + 2 > stackCapacity) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + compositeStack[stackTop] = factor; + compositeStack[stackTop + 1] = quotient; + stackTop += 2; + } + + return true; + } + + private static bool PeelSmallPrimesKernel( + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ref int factorCount, + out ulong remaining) + { + ulong remainingLocal = value; + + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > limit) + { + break; + } + + ulong square = squares[i]; + if (square != 0UL && square > remainingLocal) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remainingLocal % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); + + if (!TryAppendFactorKernel(factors, exponents, ref factorCount, primeValue, exponent)) + { + remaining = value; + return false; + } + } + + remaining = remainingLocal; + return true; + } + + private static bool TryAppendFactorKernel( + ArrayView1D factors, + ArrayView1D exponents, + ref int count, + ulong prime, + int exponent) + { + if (prime <= 1UL || exponent <= 0) + { + return true; + } + + for (int i = 0; i < count; i++) + { + if (factors[i] == prime) + { + exponents[i] += exponent; + return true; + } + } + + int capacity = (int)factors.Length; + if (count >= capacity) + { + return false; + } + + factors[count] = prime; + exponents[count] = exponent; + count++; + return true; + } + + private static ulong MulModKernel(ulong left, ulong right, ulong modulus) + { + GpuUInt128 product = new GpuUInt128(left); + return product.MulMod(right, modulus); + } + + private static ulong PollardRhoKernel(ulong value) + { + if ((value & 1UL) == 0UL) + { + return 2UL; + } + + ulong c = 1UL; + while (true) + { + ulong x = 2UL; + ulong y = 2UL; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomialKernel(x, c, value); + y = AdvancePolynomialKernel(y, c, value); + y = AdvancePolynomialKernel(y, c, value); + + ulong diff = x > y ? x - y : y - x; + d = BinaryGcdKernel(diff, value); + } + + if (d == value) + { + c++; + if (c == 0UL) + { + c = 1UL; + } + + continue; + } + + return d; + } + } + + private static ulong AdvancePolynomialKernel(ulong x, ulong c, ulong modulus) + { + ulong squared = MulModKernel(x, x, modulus); + GpuUInt128 accumulator = new GpuUInt128(squared); + accumulator.AddMod(c, modulus); + return accumulator.Low; + } + + private static ulong BinaryGcdKernel(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + ulong aLocal = a >> BitOperations.TrailingZeroCount(a); + ulong bLocal = b; + + while (true) + { + bLocal >>= BitOperations.TrailingZeroCount(bLocal); + if (aLocal > bLocal) + { + ulong temp = aLocal; + aLocal = bLocal; + bLocal = temp; + } + + bLocal -= aLocal; + if (bLocal == 0UL) + { + return aLocal << shift; + } + } + } + + private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) + { + return ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent) == 1UL; + } + + private static ulong CalculateByDoublingKernel(ulong prime) + { + ulong order = 1UL; + ulong value = 2UL % prime; + + while (value != 1UL) + { + value <<= 1; + if (value >= prime) + { + value -= prime; + } + + order++; + } + + return order; + } /// This kernel always sets the result of the corresponding element. Callers don't need to clear the output buffers. @@ -1167,7 +1206,7 @@ internal static partial class PrimeOrderGpuHeuristics remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent); } - private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) + private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) { var pool = _pow2ModWideKernel ??= []; if (pool.TryGetValue(accelerator, out var cached)) @@ -1182,194 +1221,194 @@ internal static partial class PrimeOrderGpuHeuristics cached = kernel.CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); pool[accelerator] = cached; return cached; - } + } /// This kernel always sets the result of the corresponding element. Callers don't need to clear the output buffers. - internal static void Pow2ModKernelWide(Index1D index, ArrayView1D exponents, GpuUInt128 modulus, ArrayView1D remainders) - { - GpuUInt128 exponent = exponents[index]; - remainders[index] = Pow2ModKernelCore(exponent, modulus); - } + internal static void Pow2ModKernelWide(Index1D index, ArrayView1D exponents, GpuUInt128 modulus, ArrayView1D remainders) + { + GpuUInt128 exponent = exponents[index]; + remainders[index] = Pow2ModKernelCore(exponent, modulus); + } - private static GpuUInt128[] InitializeOddPowersTable(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) - { - GpuUInt128[] result = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount]; - result[0] = baseValue; - if (oddPowerCount == 1) - { - return result; - } + private static GpuUInt128[] InitializeOddPowersTable(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + GpuUInt128[] result = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount]; + result[0] = baseValue; + if (oddPowerCount == 1) + { + return result; + } - // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. - baseValue.MulMod(baseValue, modulus); + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); - // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. - GpuUInt128 current = baseValue; + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; - // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. - for (int i = 1; i < oddPowerCount; i++) - { - current.MulMod(baseValue, modulus); - result[i] = current; - } + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + for (int i = 1; i < oddPowerCount; i++) + { + current.MulMod(baseValue, modulus); + result[i] = current; + } - return result; - } + return result; + } - internal static GpuUInt128 Pow2ModKernelCore(GpuUInt128 exponent, GpuUInt128 modulus) - { - // This should never happen in production code. - // if (modulus == GpuUInt128.One) - // { - // return GpuUInt128.Zero; - // } + internal static GpuUInt128 Pow2ModKernelCore(GpuUInt128 exponent, GpuUInt128 modulus) + { + // This should never happen in production code. + // if (modulus == GpuUInt128.One) + // { + // return GpuUInt128.Zero; + // } - // This should never happen in production code. - // if (exponent.IsZero) - // { - // return GpuUInt128.One; - // } + // This should never happen in production code. + // if (exponent.IsZero) + // { + // return GpuUInt128.One; + // } - GpuUInt128 baseValue = UInt128Numbers.TwoGpu; + GpuUInt128 baseValue = UInt128Numbers.TwoGpu; - // This should never happen in production code - 2 should never be greater or equal modulus. - // if (baseValue.CompareTo(modulus) >= 0) - // { - // baseValue.Sub(modulus); - // } + // This should never happen in production code - 2 should never be greater or equal modulus. + // if (baseValue.CompareTo(modulus) >= 0) + // { + // baseValue.Sub(modulus); + // } - if (ShouldUseSingleBit(exponent)) - { - return Pow2MontgomeryModSingleBit(exponent, modulus, baseValue); - } + if (ShouldUseSingleBit(exponent)) + { + return Pow2MontgomeryModSingleBit(exponent, modulus, baseValue); + } - int index = exponent.GetBitLength(); - int windowSize = GetWindowSize(index); - index--; + int index = exponent.GetBitLength(); + int windowSize = GetWindowSize(index); + index--; - int oddPowerCount = 1 << (windowSize - 1); - GpuUInt128[] oddPowers = InitializeOddPowersTable(baseValue, modulus, oddPowerCount); - GpuUInt128 result = GpuUInt128.One; + int oddPowerCount = 1 << (windowSize - 1); + GpuUInt128[] oddPowers = InitializeOddPowersTable(baseValue, modulus, oddPowerCount); + GpuUInt128 result = GpuUInt128.One; - int windowStart; - ulong windowValue; - while (index >= 0) - { - if (!IsBitSet(exponent, index)) - { - result.MulMod(result, modulus); - index--; - continue; - } + int windowStart; + ulong windowValue; + while (index >= 0) + { + if (!IsBitSet(exponent, index)) + { + result.MulMod(result, modulus); + index--; + continue; + } - windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } + windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } - while (!IsBitSet(exponent, windowStart)) - { - windowStart++; - } + while (!IsBitSet(exponent, windowStart)) + { + windowStart++; + } - // We're reusing oddPowerCount as windowBitCount here to lower registry pressure & avoid additional allocation - oddPowerCount = index - windowStart + 1; - index = windowStart - 1; - windowValue = ExtractWindowValue(exponent, windowStart, oddPowerCount); + // We're reusing oddPowerCount as windowBitCount here to lower registry pressure & avoid additional allocation + oddPowerCount = index - windowStart + 1; + index = windowStart - 1; + windowValue = ExtractWindowValue(exponent, windowStart, oddPowerCount); - // We're reusing windowStart as square here to lower registry pressure & avoid additional allocation - for (windowStart = 0; windowStart < oddPowerCount; windowStart++) - { - result.MulMod(result, modulus); - } + // We're reusing windowStart as square here to lower registry pressure & avoid additional allocation + for (windowStart = 0; windowStart < oddPowerCount; windowStart++) + { + result.MulMod(result, modulus); + } - // We're reusing windowStart as tableIndex here to lower registry pressure & avoid additional allocation - windowStart = (int)((windowValue - 1UL) >> 1); - result.MulMod(oddPowers[windowStart], modulus); - } + // We're reusing windowStart as tableIndex here to lower registry pressure & avoid additional allocation + windowStart = (int)((windowValue - 1UL) >> 1); + result.MulMod(oddPowers[windowStart], modulus); + } - return result; - } + return result; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool ShouldUseSingleBit(GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool ShouldUseSingleBit(GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; - private static int GetWindowSize(int bitLength) - { - if (bitLength <= Pow2WindowSizeBits) - { - return Math.Max(bitLength, 1); - } + private static int GetWindowSize(int bitLength) + { + if (bitLength <= Pow2WindowSizeBits) + { + return Math.Max(bitLength, 1); + } - if (bitLength <= 23) - { - return 4; - } + if (bitLength <= 23) + { + return 4; + } - if (bitLength <= 79) - { - return 5; - } + if (bitLength <= 79) + { + return 5; + } - if (bitLength <= 239) - { - return 6; - } + if (bitLength <= 239) + { + return 6; + } - if (bitLength <= 671) - { - return 7; - } + if (bitLength <= 671) + { + return 7; + } - return Pow2WindowSizeBits; - } + return Pow2WindowSizeBits; + } - private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) - { - GpuUInt128 result = GpuUInt128.One; + private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) + { + GpuUInt128 result = GpuUInt128.One; - while (!exponent.IsZero) - { - if ((exponent.Low & 1UL) != 0UL) - { - result.MulMod(baseValue, modulus); - } + while (!exponent.IsZero) + { + if ((exponent.Low & 1UL) != 0UL) + { + result.MulMod(baseValue, modulus); + } - exponent.ShiftRight(1); - if (exponent.IsZero) - { - break; - } + exponent.ShiftRight(1); + if (exponent.IsZero) + { + break; + } - // Reusing baseValue to store the squared base for the next iteration. - baseValue.MulMod(baseValue, modulus); - } + // Reusing baseValue to store the squared base for the next iteration. + baseValue.MulMod(baseValue, modulus); + } - return result; - } + return result; + } - private static bool IsBitSet(GpuUInt128 value, int bitIndex) - { - if (bitIndex >= 64) - { - return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; - } + private static bool IsBitSet(GpuUInt128 value, int bitIndex) + { + if (bitIndex >= 64) + { + return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; + } - return ((value.Low >> bitIndex) & 1UL) != 0UL; - } + return ((value.Low >> bitIndex) & 1UL) != 0UL; + } - private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) - { - if (windowStart != 0) - { - exponent.ShiftRight(windowStart); - ulong mask = (1UL << windowBitCount) - 1UL; - return exponent.Low & mask; - } + private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) + { + if (windowStart != 0) + { + exponent.ShiftRight(windowStart); + ulong mask = (1UL << windowBitCount) - 1UL; + return exponent.Low & mask; + } - ulong directMask = (1UL << windowBitCount) - 1UL; - return exponent.Low & directMask; - } + ulong directMask = (1UL << windowBitCount) - 1UL; + return exponent.Low & directMask; + } } diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index bcd9bda..c01c62e 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -293,7 +293,7 @@ public static class MersenneNumberDivisorByDivisorTester tasks[taskIndex++] = Task.Factory.StartNew( () => { - // startGate.Wait(); + startGate.Wait(); Console.WriteLine($"Task started for range {rangeStart}"); for (int index = rangeStart; index < rangeEnd; index++) diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index a8bf7db..9c89384 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -5,7 +5,8 @@ namespace PerfectNumbers.Core; public static class PerfectNumberConstants { public const ulong BiggestKnownEvenPerfectP = 136279841UL; - public const int ConsoleInterval = 100_000; + public const int ConsoleInterval = 100_000; + public const int DefaultFactorsBuffer = 32; public const ulong ExtraDivisorCycleSearchLimit = 64UL; public const int MaxQForDivisorCycles = 4_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 23c298a..cec27ef 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -64,7 +64,7 @@ internal static partial class PrimeOrderCalculator ulong result; if (phiFactors.Factors is null) { - result = CalculateByFactorizationCpu(gpu, prime, divisorData); + result = CalculateByFactorizationGpu(gpu, prime, divisorData); phiFactors.Dispose(); return result; } @@ -103,7 +103,7 @@ internal static partial class PrimeOrderCalculator { orderFactors?.Dispose(); GpuPrimeWorkLimiter.Release(); - return CalculateByFactorizationCpu(gpu, prime, divisorData); + return CalculateByFactorizationGpu(gpu, prime, divisorData); } if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) @@ -1658,17 +1658,17 @@ internal static partial class PrimeOrderCalculator { ulong phi = prime - 1UL; Dictionary counts = new(capacity: 8); - FactorCompletely(phi, counts); + FactorCompletelyCpu(phi, counts); if (counts.Count == 0) { return phi; } - List> entries = [.. counts]; - entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); + KeyValuePair[] entries = [.. counts]; + Array.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); ulong order = phi; - int entryCount = entries.Count; + int entryCount = entries.Length; GpuPrimeWorkLimiter.Acquire(); for (int i = 0; i < entryCount; i++) @@ -1677,12 +1677,14 @@ internal static partial class PrimeOrderCalculator int exponent = entries[i].Value; for (int iteration = 0; iteration < exponent; iteration++) { - if ((order % primeFactor) != 0UL) + ulong remainder = order % primeFactor; + if (remainder != 0UL) { break; } - ulong candidate = order / primeFactor; + ulong candidate = order - primeFactor * remainder; + // ulong candidate = order / primeFactor; // if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL) if (candidate.Pow2MontgomeryModWindowedConvertGpu(gpu, divisorData) == 1UL) @@ -1700,7 +1702,7 @@ internal static partial class PrimeOrderCalculator return order; } - private static void FactorCompletely(ulong value, Dictionary counts) + private static void FactorCompletelyCpu(ulong value, Dictionary counts) { FactorCompletely(value, counts, knownComposite: false); } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 24b066a..45fef63 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,8 +1,10 @@ using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; -using ILGPU.Runtime.Cuda; +using ILGPU.Util; +using Open.Collections; using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; @@ -12,6 +14,47 @@ internal static partial class PrimeOrderCalculator private const int GpuSmallPrimeFactorSlots = 64; + private static ulong CalculateByFactorizationGpu(Pow2MontgomeryAccelerator gpu, ulong prime, in MontgomeryDivisorData divisorData) + { + ulong phi = prime - 1UL; + Dictionary counts = gpu.ToTestOnHost; + counts.Clear(); + + FactorCompletelyCpu(phi, counts); + if (counts.Count == 0) + { + return phi; + } + + int entryCount = counts.Count; + Span> entries = stackalloc KeyValuePair[entryCount]; + + int stored = 0; + foreach (var entry in counts) + { + entries[stored++] = entry; + } + + entries.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); + + AcceleratorStream stream = gpu.Stream!; + gpu.EnsureCapacity(entryCount); + gpu.ToTestOnDevice.View.CopyFromCPU(stream, entries); + + ulong order = phi; + var kernel = gpu.CheckFactorsKernel; + stream.Synchronize(); + + GpuPrimeWorkLimiter.Acquire(); + + kernel.Launch(stream, 1, entryCount, phi, gpu.ToTestOnDevice.View, divisorData, gpu.Output.View); + gpu.Output.View.CopyToCPU(stream, ref order, 1); + stream.Synchronize(); + + GpuPrimeWorkLimiter.Release(); + return order; + } + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out int factorCount, out ulong remaining) { var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 9799a4f..773b27d 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -99,7 +99,7 @@ public sealed class PrimeTester var gpu = PrimeTesterGpuContextPool.Rent(1); var accelerator = gpu.Accelerator; var stream = gpu.Stream; - var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.GetKernel(accelerator); + var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.GetSmallPrimeSieveKernel(accelerator); ulong value = n; // Span flag = stackalloc byte[1]; @@ -167,7 +167,7 @@ public sealed class PrimeTester var gpu = PrimeTesterGpuContextPool.Rent(GpuBatchSize); var accelerator = gpu.Accelerator; var stream = gpu.Stream; - var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.GetKernel(accelerator); + var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.GetSmallPrimeSieveKernel(accelerator); int totalLength = values.Length; int batchSize = GpuBatchSize; @@ -324,11 +324,11 @@ public sealed class PrimeTester } [ThreadStatic] - public static Dictionary, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>? Kernel; + public static Dictionary, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>? SmallPrimeSieveKernel; - internal static Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> GetKernel(Accelerator accelerator) + internal static Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> GetSmallPrimeSieveKernel(Accelerator accelerator) { - var pool = Kernel ??= []; + var pool = SmallPrimeSieveKernel ??= []; if (pool.TryGetValue(accelerator, out var cached)) { return cached;