From 26cc57402035e7fdbde8216d9dfed27065b25201 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 5 Oct 2025 06:26:09 +0200 Subject: [PATCH] Treat Pollard Rho success as the GPU primality check --- PerfectNumbers.Core/DivisorCycleCache.cs | 30 + .../Gpu/PrimeOrderGpuHeuristics.cs | 1521 ++++++++++++++++- PerfectNumbers.Core/MersenneDivisorCycles.cs | 30 +- .../PrimeOrderCalculator.Wide.cs | 135 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 654 +++++-- todo.txt | 1 - tools/GpuNumericTypeProbe.csx | 30 +- 7 files changed, 2188 insertions(+), 213 deletions(-) diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 8e0b61d..0da555a 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -262,9 +262,39 @@ public sealed class DivisorCycleCache } ulong divisor = divisors[index]; + + if ((divisor & (divisor - 1UL)) == 0UL) + { + cycles[index] = 1UL; + status[index] = ByteOne; + pow[index] = 1UL; + order[index] = 1UL; + return; + } + ulong currentPow = pow[index]; ulong currentOrder = order[index]; + if (divisor <= 3UL) + { + while (currentPow != 1UL) + { + currentPow += currentPow; + if (currentPow >= divisor) + { + currentPow -= divisor; + } + + currentOrder++; + } + + cycles[index] = currentOrder; + status[index] = ByteOne; + pow[index] = currentPow; + order[index] = currentOrder; + return; + } + do { currentPow += currentPow; diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index e3da840..5bd2eab 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -16,18 +16,114 @@ internal enum GpuPow2ModStatus Unavailable, } -internal static class PrimeOrderGpuHeuristics +internal static partial class PrimeOrderGpuHeuristics { private static readonly ConcurrentDictionary OverflowedPrimes = new(); private static readonly ConcurrentDictionary OverflowedPrimesWide = new(); private static readonly ConcurrentDictionary, MontgomeryDivisorData, ArrayView1D>> Pow2ModKernelCache = new(); private static readonly ConcurrentDictionary, GpuUInt128, ArrayView1D>> Pow2ModKernelWideCache = new(); + private static readonly ConcurrentDictionary, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>> PartialFactorKernelCache = new(); + private static readonly ConcurrentDictionary OrderKernelCache = new(); + private static readonly ConcurrentDictionary SmallPrimeDeviceCaches = new(); + + private readonly struct OrderKernelConfig + { + public OrderKernelConfig(ulong previousOrder, byte hasPreviousOrder, uint smallFactorLimit, int maxPowChecks, int mode) + { + PreviousOrder = previousOrder; + HasPreviousOrder = hasPreviousOrder; + SmallFactorLimit = smallFactorLimit; + MaxPowChecks = maxPowChecks; + Mode = mode; + } + + public ulong PreviousOrder { get; } + + public byte HasPreviousOrder { get; } + + public uint SmallFactorLimit { get; } + + public int MaxPowChecks { get; } + + public int Mode { get; } + } + + private readonly struct OrderKernelBuffers + { + public OrderKernelBuffers( + ArrayView1D phiFactors, + ArrayView1D phiExponents, + ArrayView1D workFactors, + ArrayView1D workExponents, + ArrayView1D candidates, + ArrayView1D stackIndex, + ArrayView1D stackExponent, + ArrayView1D stackProduct, + ArrayView1D result, + ArrayView1D status) + { + PhiFactors = phiFactors; + PhiExponents = phiExponents; + WorkFactors = workFactors; + WorkExponents = workExponents; + Candidates = candidates; + StackIndex = stackIndex; + StackExponent = stackExponent; + StackProduct = stackProduct; + Result = result; + Status = status; + } + + public ArrayView1D PhiFactors { get; } + + public ArrayView1D PhiExponents { get; } + + public ArrayView1D WorkFactors { get; } + + public ArrayView1D WorkExponents { get; } + + public ArrayView1D Candidates { get; } + + public ArrayView1D StackIndex { get; } + + public ArrayView1D StackExponent { get; } + + public ArrayView1D StackProduct { get; } + + public ArrayView1D Result { get; } + + public ArrayView1D Status { get; } + } + + private delegate void OrderKernelLauncher( + AcceleratorStream stream, + Index1D extent, + ulong prime, + OrderKernelConfig config, + MontgomeryDivisorData divisor, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + OrderKernelBuffers buffers); + + private sealed class SmallPrimeDeviceCache + { + public MemoryBuffer1D? Primes; + public MemoryBuffer1D? Squares; + public int Count; + } + private const int WideStackThreshold = 8; private static PrimeOrderGpuCapability s_capability = PrimeOrderGpuCapability.Default; private const int Pow2WindowSizeBits = 8; private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); private const ulong Pow2WindowFallbackThreshold = 32UL; + private const int HeuristicCandidateLimit = 512; + private const int HeuristicStackCapacity = 256; + + private const int GpuSmallPrimeFactorSlots = 64; + [InlineArray(Pow2WindowOddPowerCount)] private struct Pow2OddPowerTable @@ -48,6 +144,223 @@ internal static class PrimeOrderGpuHeuristics s_capability = PrimeOrderGpuCapability.Default; } + private static SmallPrimeDeviceCache GetSmallPrimeDeviceCache(Accelerator accelerator) + { + return SmallPrimeDeviceCaches.GetOrAdd(accelerator, static acc => + { + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + var primeBuffer = acc.Allocate1D(primes.Length); + primeBuffer.View.CopyFromCPU(primes); + var squareBuffer = acc.Allocate1D(squares.Length); + squareBuffer.View.CopyFromCPU(squares); + return new SmallPrimeDeviceCache + { + Primes = primeBuffer, + Squares = squareBuffer, + Count = primes.Length, + }; + }); + } + + private static Action, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> GetPartialFactorKernel(Accelerator accelerator) + { + return PartialFactorKernelCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(PartialFactorKernel); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + }); + } + + public static bool TryPartialFactor( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining, + out bool fullyFactored) + { + factorCount = 0; + remaining = value; + fullyFactored = false; + + if (primeTargets.Length == 0 || exponentTargets.Length == 0) + { + return false; + } + + primeTargets.Clear(); + exponentTargets.Clear(); + + if (!TryLaunchPartialFactorKernel( + value, + limit, + primeTargets, + exponentTargets, + out int extracted, + out ulong leftover, + out bool kernelFullyFactored)) + { + return false; + } + + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + if (extracted > capacity) + { + factorCount = 0; + remaining = value; + fullyFactored = false; + return false; + } + + factorCount = extracted; + remaining = leftover; + fullyFactored = kernelFullyFactored && leftover == 1UL; + return true; + } + + private static bool TryLaunchPartialFactorKernel( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining, + out bool fullyFactored) + { + factorCount = 0; + remaining = value; + fullyFactored = false; + + try + { + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + try + { + using var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + var kernel = GetPartialFactorKernel(accelerator); + SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); + + using var factorBuffer = accelerator.Allocate1D(primeTargets.Length); + using var exponentBuffer = accelerator.Allocate1D(exponentTargets.Length); + using var countBuffer = accelerator.Allocate1D(1); + using var remainingBuffer = accelerator.Allocate1D(1); + using var fullyFactoredBuffer = accelerator.Allocate1D(1); + + factorBuffer.MemSetToZero(); + exponentBuffer.MemSetToZero(); + countBuffer.MemSetToZero(); + remainingBuffer.MemSetToZero(); + fullyFactoredBuffer.MemSetToZero(); + + kernel( + stream, + 1, + cache.Primes!.View, + cache.Squares!.View, + cache.Count, + primeTargets.Length, + value, + limit, + factorBuffer.View, + exponentBuffer.View, + countBuffer.View, + remainingBuffer.View, + fullyFactoredBuffer.View); + + stream.Synchronize(); + + countBuffer.View.CopyToCPU(ref factorCount, 1); + factorCount = Math.Min(factorCount, primeTargets.Length); + factorBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(primeTargets), primeTargets.Length); + exponentBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentTargets), exponentTargets.Length); + remainingBuffer.View.CopyToCPU(ref remaining, 1); + + byte fullyFactoredFlag = 0; + fullyFactoredBuffer.View.CopyToCPU(ref fullyFactoredFlag, 1); + fullyFactored = fullyFactoredFlag != 0; + + return true; + } + finally + { + lease.Dispose(); + } + } + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + { + factorCount = 0; + remaining = value; + fullyFactored = false; + return false; + } + } + + private 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; + + 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 primeValue = primeCandidate; + if (primeValue == 0UL || (remainingLocal % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); + + factorsOut[count] = primeValue; + exponentsOut[count] = exponent; + count++; + } + + countOut[0] = count; + remainingOut[0] = remainingLocal; + fullyFactoredOut[0] = remainingLocal == 1UL ? (byte)1 : (byte)0; + } public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder) { Span exponents = stackalloc ulong[1]; @@ -118,6 +431,1212 @@ internal static class PrimeOrderGpuHeuristics } public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan exponents, UInt128 prime, Span remainders) + { + return TryPow2ModBatchInternal(exponents, prime, remainders); + } + + internal static bool TryCalculateOrder( + ulong prime, + ulong? previousOrder, + PrimeOrderCalculator.PrimeOrderSearchConfig config, + in MontgomeryDivisorData divisorData, + out PrimeOrderCalculator.PrimeOrderResult result) + { + result = default; + + try + { + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + try + { + using var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + var kernel = GetOrderKernel(accelerator); + SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); + + using var phiFactorBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + using var phiExponentBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + using var workFactorBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + using var workExponentBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + using var candidateBuffer = accelerator.Allocate1D(HeuristicCandidateLimit); + using var stackIndexBuffer = accelerator.Allocate1D(HeuristicStackCapacity); + using var stackExponentBuffer = accelerator.Allocate1D(HeuristicStackCapacity); + using var stackProductBuffer = accelerator.Allocate1D(HeuristicStackCapacity); + using var resultBuffer = accelerator.Allocate1D(1); + using var statusBuffer = accelerator.Allocate1D(1); + + phiFactorBuffer.MemSetToZero(); + phiExponentBuffer.MemSetToZero(); + workFactorBuffer.MemSetToZero(); + workExponentBuffer.MemSetToZero(); + candidateBuffer.MemSetToZero(); + stackIndexBuffer.MemSetToZero(); + stackExponentBuffer.MemSetToZero(); + stackProductBuffer.MemSetToZero(); + resultBuffer.MemSetToZero(); + statusBuffer.MemSetToZero(); + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + byte hasPrevious = previousOrder.HasValue ? (byte)1 : (byte)0; + ulong previousValue = previousOrder ?? 0UL; + + var kernelConfig = new OrderKernelConfig(previousValue, hasPrevious, limit, config.MaxPowChecks, (int)config.Mode); + var buffers = new OrderKernelBuffers( + phiFactorBuffer.View, + phiExponentBuffer.View, + workFactorBuffer.View, + workExponentBuffer.View, + candidateBuffer.View, + stackIndexBuffer.View, + stackExponentBuffer.View, + stackProductBuffer.View, + resultBuffer.View, + statusBuffer.View); + + kernel( + stream, + 1, + prime, + kernelConfig, + divisorData, + cache.Primes!.View, + cache.Squares!.View, + cache.Count, + buffers); + + stream.Synchronize(); + + byte status = 0; + statusBuffer.View.CopyToCPU(ref status, 1); + + PrimeOrderKernelStatus kernelStatus = (PrimeOrderKernelStatus)status; + if (kernelStatus == PrimeOrderKernelStatus.Fallback) + { + return false; + } + + if (kernelStatus == PrimeOrderKernelStatus.PollardOverflow) + { + throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity."); + } + + ulong order = 0UL; + resultBuffer.View.CopyToCPU(ref order, 1); + + PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found + ? PrimeOrderCalculator.PrimeOrderStatus.Found + : PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved; + + result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order); + return true; + } + finally + { + lease.Dispose(); + } + } + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + { + result = default; + return false; + } + } + + private static OrderKernelLauncher GetOrderKernel(Accelerator accelerator) + { + return OrderKernelCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(CalculateOrderKernel); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate(); + }); + } + + private 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; + } + + public int Primary { get; } + + public long Secondary { get; } + + public long Tertiary { get; } + + public int CompareTo(CandidateKey other) + { + if (Primary != other.Primary) + { + return Primary.CompareTo(other.Primary); + } + + if (Secondary != other.Secondary) + { + return Secondary.CompareTo(other.Secondary); + } + + return Tertiary.CompareTo(other.Tertiary); + } + } + + private static void CalculateOrderKernel( + Index1D index, + ulong prime, + OrderKernelConfig config, + MontgomeryDivisorData divisor, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + OrderKernelBuffers buffers) + { + if (index != 0) + { + return; + } + + 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; + + 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; + + if (prime <= 3UL) + { + ulong orderValue = prime == 3UL ? 2UL : 1UL; + resultOut[0] = orderValue; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + 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 (Pow2EqualsOneKernel(reduced, prime, divisor)) + { + 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, prime, 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, prime, 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, prime, 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, prime, 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, prime, 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, prime, 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, ulong prime, in MontgomeryDivisorData divisor) + { + return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 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; + } + + private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan exponents, UInt128 prime, Span remainders) { if (exponents.Length == 0) { diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 1de140d..2e01d4f 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -143,7 +143,8 @@ public class MersenneDivisorCycles PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate( divisor, previousOrder: null, - PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); if (wideResult.Order != UInt128.Zero) { return wideResult.Order; @@ -192,7 +193,8 @@ public class MersenneDivisorCycles PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( divisor, previousOrder: null, - PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); if (orderResult.Order != 0UL) { cycleLength = orderResult.Order; @@ -818,15 +820,24 @@ public class MersenneDivisorCycles public static ulong CalculateCycleLength(ulong divisor) { - // For divisor = 2^k, cycle is 1 + if (TryCalculateCycleLengthHeuristic(divisor, out ulong cycleLength)) + return cycleLength; + + return CalculateCycleLengthFallback(divisor); + } + + internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, out ulong cycleLength) + { if ((divisor & (divisor - 1UL)) == 0UL) { - return 1UL; + cycleLength = 1UL; + return true; } if (divisor <= 3UL) { - return CalculateCycleLengthFallback(divisor); + cycleLength = CalculateCycleLengthFallback(divisor); + return true; } if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) @@ -834,14 +845,17 @@ public class MersenneDivisorCycles PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( divisor, previousOrder: null, - PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); if (orderResult.Order != 0UL) { - return orderResult.Order; + cycleLength = orderResult.Order; + return true; } } - return CalculateCycleLengthFallback(divisor); + cycleLength = 0UL; + return false; } private static ulong CalculateCycleLengthFallback(ulong divisor) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs index 9e85e7b..3d618ba 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs @@ -23,7 +23,15 @@ internal static partial class PrimeOrderCalculator } public static PrimeOrderResultWide Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config) + => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu); + + public static PrimeOrderResultWide Calculate( + UInt128 prime, + UInt128? previousOrder, + PrimeOrderSearchConfig config, + PrimeOrderHeuristicDevice device) { + using var scope = UsePow2Mode(device); if (prime <= ulong.MaxValue) { ulong? previous = null; @@ -40,7 +48,7 @@ internal static partial class PrimeOrderCalculator } } - PrimeOrderResult result = Calculate((ulong)prime, previous, config); + PrimeOrderResult result = Calculate((ulong)prime, previous, config, device); return new PrimeOrderResultWide(result.Status, result.Order); } @@ -60,12 +68,23 @@ internal static partial class PrimeOrderCalculator } UInt128 phi = prime - UInt128.One; + using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice); PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); if (phiFactors.Factors is null) { return FinishStrictlyWide(prime, config.Mode); } + return RunHeuristicPipelineWide(prime, previousOrder, config, phi, phiFactors); + } + + private static PrimeOrderResultWide RunHeuristicPipelineWide( + UInt128 prime, + UInt128? previousOrder, + PrimeOrderSearchConfig config, + UInt128 phi, + PartialFactorResult128 phiFactors) + { if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, phiFactors)) { return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi); @@ -553,41 +572,29 @@ internal static partial class PrimeOrderCalculator } Dictionary counts = new(capacity: 8); - UInt128 remaining = value; - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + UInt128 remaining; - for (int i = 0; i < primeCount; i++) + if (IsGpuPow2Allowed && value <= ulong.MaxValue) { - uint primeCandidate = primes[i]; - if (primeCandidate > limit) + Dictionary narrowCounts = new(capacity: 8); + if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining)) { - break; - } + foreach (KeyValuePair entry in narrowCounts) + { + counts[(UInt128)entry.Key] = entry.Value; + } - UInt128 primeValue = primeCandidate; - UInt128 primeSquare = (UInt128)squares[i]; - if (primeSquare > remaining) + remaining = (UInt128)narrowRemaining; + } + else { - break; + remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); } - - if ((remaining % primeValue) != UInt128.Zero) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == UInt128.Zero); - - counts[primeValue] = exponent; + } + else + { + remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); } List pending = new(); @@ -682,6 +689,46 @@ internal static partial class PrimeOrderCalculator return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index); } + private static UInt128 PopulateSmallPrimeFactorsCpuWide(UInt128 value, uint limit, Dictionary counts) + { + UInt128 remaining = value; + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + + for (int i = 0; i < primeCount; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > limit) + { + break; + } + + UInt128 primeSquare = (UInt128)squares[i]; + if (primeSquare > remaining) + { + break; + } + + UInt128 primeValue = primeCandidate; + if ((remaining % primeValue) != UInt128.Zero) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == UInt128.Zero); + + counts[primeValue] = exponent; + } + + return remaining; + } private static bool TryPollardRhoWide(UInt128 n, Stopwatch stopwatch, long budgetTicks, out UInt128 factor) { factor = UInt128.Zero; @@ -955,21 +1002,25 @@ internal static partial class PrimeOrderCalculator { return UInt128.Zero; } - if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue) - { - ulong prime64 = (ulong)modulus; - ulong exponent64 = (ulong)exponent; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder); - if (status == GpuPow2ModStatus.Success) - { - return remainder; - } - } - GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, modulus, out UInt128 wideRemainder); - if (wideStatus == GpuPow2ModStatus.Success) + if (IsGpuPow2Allowed) { - return wideRemainder; + if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue) + { + ulong prime64 = (ulong)modulus; + ulong exponent64 = (ulong)exponent; + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder); + if (status == GpuPow2ModStatus.Success) + { + return remainder; + } + } + + GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, modulus, out UInt128 wideRemainder); + if (wideStatus == GpuPow2ModStatus.Success) + { + return wideRemainder; + } } return exponent.Pow2MontgomeryModWindowed(modulus); diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 240d4ed..7541403 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -1,3 +1,4 @@ +using System; using System.Buffers; using System.Buffers.Binary; using System.Collections.Generic; @@ -61,7 +62,104 @@ internal static partial class PrimeOrderCalculator public static PrimeOrderSearchConfig StrictDefault => new(1_000_000, 0, 0, PrimeOrderMode.Strict); } + internal enum PrimeOrderHeuristicDevice + { + Cpu, + Gpu, + } + + [ThreadStatic] + private static bool s_pow2ModeInitialized; + + [ThreadStatic] + private static bool s_allowGpuPow2; + + private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; + + + [ThreadStatic] + private static PrimeOrderHeuristicDevice s_deviceMode; + + [ThreadStatic] + private static bool s_debugLoggingEnabled; + + private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; + private readonly struct Pow2ModeScope : IDisposable + { + private readonly bool _previousInitialized; + private readonly bool _previousAllow; + private readonly PrimeOrderHeuristicDevice _previousDevice; + + public Pow2ModeScope(PrimeOrderHeuristicDevice device) + { + _previousInitialized = s_pow2ModeInitialized; + _previousAllow = s_allowGpuPow2; + _previousDevice = s_deviceMode; + s_pow2ModeInitialized = true; + s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; + s_deviceMode = device; + } + + public void Dispose() + { + s_allowGpuPow2 = _previousAllow; + s_pow2ModeInitialized = _previousInitialized; + s_deviceMode = _previousDevice; + } + } + + private readonly struct DebugLoggingScope : IDisposable + { + private readonly bool _previous; + + public DebugLoggingScope(bool enabled) + { + _previous = s_debugLoggingEnabled; + s_debugLoggingEnabled = enabled; + } + + public void Dispose() + { + s_debugLoggingEnabled = _previous; + } + } + + private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device); + + private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled); + + [Conditional("DEBUG")] + private static void DebugLog(string message) + { + if (s_debugLoggingEnabled) + { + Console.WriteLine(message); + } + } + + [Conditional("DEBUG")] + private static void DebugLog(Func messageFactory) + { + if (s_debugLoggingEnabled) + { + Console.WriteLine(messageFactory()); + } + } + public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config) + => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu); + + public static PrimeOrderResult Calculate( + ulong prime, + ulong? previousOrder, + PrimeOrderSearchConfig config, + PrimeOrderHeuristicDevice device) + { + using var scope = UsePow2Mode(device); + return CalculateInternal(prime, previousOrder, config); + } + + private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config) { if (prime <= 3UL) { @@ -71,36 +169,44 @@ internal static partial class PrimeOrderCalculator ulong phi = prime - 1UL; MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(prime); -#if DEBUG - Console.WriteLine("Partial factoring φ(p)"); -#endif + if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult)) + { + return gpuResult; + } + + using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice); + DebugLog("Partial factoring φ(p)"); + PartialFactorResult phiFactors = PartialFactor(phi, config); if (phiFactors.Factors is null) { -#if DEBUG - Console.WriteLine("No factors found"); -#endif + DebugLog("No factors found"); return FinishStrictly(prime, config.Mode); } -#if DEBUG - Console.WriteLine("Trying special max check"); -#endif + return RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors); + } + + private static PrimeOrderResult RunHeuristicPipeline( + ulong prime, + ulong? previousOrder, + PrimeOrderSearchConfig config, + in MontgomeryDivisorData divisorData, + ulong phi, + PartialFactorResult phiFactors) + { + DebugLog("Trying special max check"); if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData)) { return new PrimeOrderResult(PrimeOrderStatus.Found, phi); } -#if DEBUG - Console.WriteLine("Initializing starting order"); -#endif + DebugLog("Initializing starting order"); ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData); candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData); -#if DEBUG - Console.WriteLine("Trying to confirm order"); -#endif + DebugLog("Trying to confirm order"); if (TryConfirmOrder(prime, candidateOrder, divisorData, config)) { return new PrimeOrderResult(PrimeOrderStatus.Found, candidateOrder); @@ -116,9 +222,7 @@ internal static partial class PrimeOrderCalculator return new PrimeOrderResult(PrimeOrderStatus.Found, order); } -#if DEBUG - Console.WriteLine("Heuristic unresolved, finishing strictly"); -#endif + DebugLog("Heuristic unresolved, finishing strictly"); return FinishStrictly(prime, config.Mode); } @@ -222,17 +326,13 @@ internal static partial class PrimeOrderCalculator } // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`. -#if DEBUG - Console.WriteLine("Verifying a^order ≡ 1 (mod p)"); -#endif + DebugLog("Verifying a^order ≡ 1 (mod p)"); if (!Pow2EqualsOne(order, prime, divisorData)) { return false; } -#if DEBUG - Console.WriteLine("Partial factoring order"); -#endif + DebugLog("Partial factoring order"); // TODO: Do we do partial factoring of order multiple times? PartialFactorResult factorization = PartialFactor(order, config); @@ -245,16 +345,12 @@ internal static partial class PrimeOrderCalculator { if (factorization.Cofactor <= 1UL) { -#if DEBUG - Console.WriteLine("Cofactor <= 1. No factors found"); -#endif + DebugLog("Cofactor <= 1. No factors found"); return false; } // TODO: Use Open.Numerics.Primality for this final check once it's available. -#if DEBUG - Console.WriteLine("Cofactor > 1. Testing primality of cofactor"); -#endif + DebugLog("Cofactor > 1. Testing primality of cofactor"); bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); if (!isPrime) @@ -262,16 +358,12 @@ internal static partial class PrimeOrderCalculator return false; } -#if DEBUG - Console.WriteLine("Adding cofactor as prime factor"); -#endif + DebugLog("Adding cofactor as prime factor"); factorization = factorization.WithAdditionalPrime(factorization.Cofactor); } ReadOnlySpan span = factorization.Factors; -#if DEBUG - Console.WriteLine("Verifying prime-power reductions"); -#endif + DebugLog("Verifying prime-power reductions"); int length = factorization.Count; for (int i = 0; i < length; i++) { @@ -311,9 +403,7 @@ internal static partial class PrimeOrderCalculator } // TODO: Do we do partial factoring of order multiple times? -#if DEBUG - Console.WriteLine("Trying heuristic. Partial factoring order"); -#endif + DebugLog("Trying heuristic. Partial factoring order"); PartialFactorResult orderFactors = PartialFactor(order, config); if (orderFactors.Factors is null) { @@ -340,18 +430,14 @@ internal static partial class PrimeOrderCalculator int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; List candidates = new(capacity); FactorEntry[] factorArray = orderFactors.Factors!; -#if DEBUG - Console.WriteLine("Building candidates list"); -#endif + DebugLog("Building candidates list"); BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); if (candidates.Count == 0) { return false; } -#if DEBUG - Console.WriteLine("Sorting candidates"); -#endif + DebugLog("Sorting candidates"); SortCandidates(prime, previousOrder, candidates); int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; @@ -360,9 +446,7 @@ internal static partial class PrimeOrderCalculator bool allowGpuBatch = true; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); -#if DEBUG - Console.WriteLine($"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); -#endif + DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); int index = 0; const int MaxGpuBatchSize = 256; const int StackGpuBatchSize = 32; @@ -384,7 +468,7 @@ internal static partial class PrimeOrderCalculator bool gpuStackRemainders = false; GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; - if (allowGpuBatch) + if (allowGpuBatch && IsGpuPow2Allowed) { if (batchSize <= StackGpuBatchSize) { @@ -461,9 +545,7 @@ internal static partial class PrimeOrderCalculator index += batchSize; } -#if DEBUG - Console.WriteLine("No candidate confirmed"); -#endif + DebugLog("No candidate confirmed"); return false; } @@ -691,10 +773,13 @@ internal static partial class PrimeOrderCalculator private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder); - if (status == GpuPow2ModStatus.Success) + if (IsGpuPow2Allowed) { - return remainder == 1UL; + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder); + if (status == GpuPow2ModStatus.Success) + { + return remainder == 1UL; + } } return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; @@ -707,20 +792,261 @@ internal static partial class PrimeOrderCalculator return PartialFactorResult.Empty; } - Dictionary counts = new(capacity: 8); + const int FactorSlotCount = GpuSmallPrimeFactorSlots; + Span primeSlots = stackalloc ulong[FactorSlotCount]; + Span exponentSlots = stackalloc int[FactorSlotCount]; + primeSlots.Clear(); + exponentSlots.Clear(); + + int factorCount = 0; + Dictionary? counts = null; + bool useDictionary = false; + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong remaining; + + bool gpuFactored = false; + if (IsGpuHeuristicDevice) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); + } + else + { + remaining = value; + } + + if (!gpuFactored) + { + primeSlots.Clear(); + exponentSlots.Clear(); + if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) + { + useDictionary = true; + counts = new Dictionary(capacity: 8); + remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + } + } + + List pending = new(); + if (remaining > 1UL) + { + pending.Add(remaining); + } + + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + DebugLog("Processing pending composites with Pollard's Rho"); + Stopwatch stopwatch = Stopwatch.StartNew(); + long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks; + Stack stack = new(); + stack.Push(remaining); + pending.Clear(); + + while (stack.Count > 0) + { + ulong composite = stack.Pop(); + if (composite <= 1UL) + { + continue; + } + + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + continue; + } + + if (stopwatch.ElapsedTicks > budgetTicks) + { + pending.Add(composite); + continue; + } + + if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor)) + { + pending.Add(composite); + continue; + } + + ulong quotient = composite / factor; + stack.Push(factor); + stack.Push(quotient); + } + } + + ulong cofactor = 1UL; + int pendingCount = pending.Count; + for (int i = 0; i < pendingCount; i++) + { + ulong composite = pending[i]; + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + } + else if (composite > 1UL) + { + cofactor = checked(cofactor * composite); + } + } + + if (useDictionary) + { + if ((counts is null || counts.Count == 0) && cofactor == value) + { + DebugLog("cofactor is the same as value, no factors found"); + return new PartialFactorResult(null, value, false, 0); + } + + if (counts is null || counts.Count == 0) + { + return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0); + } + + FactorEntry[] factors = ArrayPool.Shared.Rent(counts.Count); + int index = 0; + DebugLog(() => $"Collecting {counts.Count} prime factors"); + 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)); + + FactorEntry[] resultArray = new FactorEntry[index]; + span.CopyTo(resultArray); + ArrayPool.Shared.Return(factors, clearArray: false); + + bool fullyFactored = cofactor == 1UL; + return new PartialFactorResult(resultArray, 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) + { + DebugLog("cofactor is the same as value, no factors found"); + return new PartialFactorResult(null, value, false, 0); + } + + return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0); + } + + FactorEntry[] array = new FactorEntry[actualCount]; + int arrayIndex = 0; + DebugLog(() => $"Collecting {actualCount} prime factors"); + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeSlots[i]; + int exponentValue = exponentSlots[i]; + if (primeValue == 0UL || exponentValue == 0) + { + continue; + } + + array[arrayIndex] = new FactorEntry(primeValue, exponentValue); + arrayIndex++; + } + + Array.Sort(array, static (a, b) => a.Value.CompareTo(b.Value)); + bool fullyFactoredArray = cofactor == 1UL; + return new PartialFactorResult(array, cofactor, fullyFactoredArray, array.Length); + } + + private const int GpuSmallPrimeFactorSlots = 64; + + private static bool TryPopulateSmallPrimeFactorsCpu( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) + { + factorCount = 0; + remaining = value; + + if (value <= 1UL) + { + return true; + } + + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + if (capacity == 0) + { + return false; + } + + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + ulong remainingLocal = value; + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remainingLocal) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remainingLocal % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); + + if (factorCount >= capacity) + { + return false; + } + + primeTargets[factorCount] = primeValue; + exponentTargets[factorCount] = exponent; + factorCount++; + } + + remaining = remainingLocal; + return true; + } + + private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) + { ulong remaining = value; uint[] primes = PrimesGenerator.SmallPrimes; ulong[] squares = PrimesGenerator.SmallPrimesPow2; int primeCount = primes.Length; - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; for (int i = 0; i < primeCount; i++) { uint primeCandidate = primes[i]; - -#if DEBUG - Console.WriteLine($"Trying small prime {primeCandidate}"); -#endif if (primeCandidate > limit) { break; @@ -748,107 +1074,37 @@ internal static partial class PrimeOrderCalculator counts[primeValue] = exponent; } - List pending = new(); - if (remaining > 1UL) - { - pending.Add(remaining); - } - - if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) - { -#if DEBUG - Console.WriteLine("Processing pending composites with Pollard's Rho"); -#endif - Stopwatch stopwatch = Stopwatch.StartNew(); - long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks; - Stack stack = new(); - stack.Push(remaining); - pending.Clear(); - - while (stack.Count > 0) - { - ulong composite = stack.Pop(); - if (composite == 1UL) - { - continue; - } - - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - // if (PrimeTester.IsPrimeInternal(composite, CancellationToken.None)) - if (isPrime) - { - AddFactor(counts, composite, 1); - continue; - } - - if (stopwatch.ElapsedTicks > budgetTicks) - { - pending.Add(composite); - continue; - } - - if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor)) - { - pending.Add(composite); - continue; - } - - ulong quotient = composite / factor; - stack.Push(factor); - stack.Push(quotient); - } - } - - ulong cofactor = 1UL; - int pendingCount = pending.Count; -#if DEBUG - Console.WriteLine($"Processing {pendingCount} pending composites with Open.Numeric.Primes"); -#endif - for (int i = 0; i < pendingCount; i++) - { - ulong composite = pending[i]; - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - // if (PrimeTester.IsPrimeInternal(composite, CancellationToken.None)) - if (isPrime) - { - AddFactor(counts, composite, 1); - } - else - { - cofactor = checked(cofactor * composite); - } - } - - if (counts.Count == 0 && cofactor == value) - { -#if DEBUG - Console.WriteLine("cofactor is the same as value, no factors found"); -#endif - return new PartialFactorResult(null, value, false, 0); - } - - FactorEntry[] factors = ArrayPool.Shared.Rent(counts.Count); - int index = 0; -#if DEBUG - Console.WriteLine($"Collecting {counts.Count} prime factors"); -#endif - 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)); - - FactorEntry[] resultArray = new FactorEntry[index]; - span.CopyTo(resultArray); - ArrayPool.Shared.Return(factors, clearArray: false); - - bool fullyFactored = cofactor == 1UL; - return new PartialFactorResult(resultArray, cofactor, fullyFactored, index); + return remaining; } + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) + { + Span primeBuffer = stackalloc ulong[GpuSmallPrimeFactorSlots]; + Span exponentBuffer = stackalloc int[GpuSmallPrimeFactorSlots]; + primeBuffer.Clear(); + exponentBuffer.Clear(); + + if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) + { + remaining = value; + return false; + } + + remaining = gpuRemaining; + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + if (primeValue == 0UL || exponent == 0) + { + continue; + } + + counts[primeValue] = exponent; + } + + return true; + } private static bool TryPollardRho(ulong n, Stopwatch stopwatch, long budgetTicks, out ulong factor) { factor = 0UL; @@ -946,6 +1202,84 @@ internal static partial class PrimeOrderCalculator } } + private static bool TryAppendFactor(Span primes, Span exponents, ref int count, ulong prime, int exponent) + { + if (exponent <= 0) + { + return true; + } + + for (int i = 0; i < count; i++) + { + if (primes[i] == prime) + { + exponents[i] += exponent; + return true; + } + } + + if (count >= primes.Length) + { + return false; + } + + primes[count] = prime; + exponents[count] = exponent; + count++; + return true; + } + + private static void AddFactorToCollector( + ref bool useDictionary, + ref Dictionary? counts, + Span primes, + Span exponents, + ref int count, + ulong prime, + int exponent) + { + if (exponent <= 0) + { + return; + } + + if (useDictionary) + { + AddFactor(counts!, prime, exponent); + return; + } + + if (TryAppendFactor(primes, exponents, ref count, prime, exponent)) + { + return; + } + + counts ??= new Dictionary(Math.Max(count, 8)); + CopyFactorsToDictionary(primes, exponents, count, counts); + count = 0; + useDictionary = true; + AddFactor(counts, prime, exponent); + } + + private static void CopyFactorsToDictionary( + ReadOnlySpan primes, + ReadOnlySpan exponents, + int count, + Dictionary target) + { + for (int i = 0; i < count; i++) + { + ulong prime = primes[i]; + int exponent = exponents[i]; + if (prime == 0UL || exponent == 0) + { + continue; + } + + target[prime] = exponent; + } + } + private static ulong CalculateByDoubling(ulong prime) { ulong order = 1UL; diff --git a/todo.txt b/todo.txt index 4acde84..19278dc 100644 --- a/todo.txt +++ b/todo.txt @@ -1,2 +1 @@ -- Mirror the complete CPU heuristic cycle algorithm on the GPU path so both implementations stay aligned. - Extend the GPU numeric type probe to cover proposed custom types such as GpuUInt256 and GpuUInt512 once their definitions are available so we can plan the wide-order heuristics on the device. diff --git a/tools/GpuNumericTypeProbe.csx b/tools/GpuNumericTypeProbe.csx index 0bfb2a6..fe3d284 100644 --- a/tools/GpuNumericTypeProbe.csx +++ b/tools/GpuNumericTypeProbe.csx @@ -1,5 +1,6 @@ #r "nuget: ILGPU, 1.5.3" #r "nuget: PeterO.Numbers, 1.8.1" +#r "nuget: Open.Numeric.Primes, 4.0.4" using System; using System.Numerics; @@ -7,6 +8,7 @@ using ILGPU; using ILGPU.Runtime; using ILGPU.Runtime.CPU; using PeterO.Numbers; +using Open.Numeric.Primes; var context = Context.CreateDefault(); var accelerator = context.CreateCPUAccelerator(0); @@ -20,7 +22,8 @@ try [ ProbeBigInteger(accelerator), ProbeEInteger(accelerator), - ProbeERational(accelerator) + ProbeERational(accelerator), + ProbeOpenNumericPrimesIsPrime(accelerator) ]; foreach (ProbeResult result in results) @@ -112,6 +115,31 @@ static ProbeResult ProbeERational(Accelerator accelerator) } } +static ProbeResult ProbeOpenNumericPrimesIsPrime(Accelerator accelerator) +{ + try + { + var kernel = accelerator.LoadAutoGroupedStreamKernel>(OpenNumericPrimesIsPrimeKernel); + using var buffer = accelerator.Allocate1D(4); + kernel(new Index1D((int)buffer.Length), buffer.View); + accelerator.Synchronize(); + return ProbeResult.Success( + "Open.Numeric.Primes.Prime.Numbers.IsPrime", + "IsPrime compiles and runs on the CPU accelerator kernel."); + } + catch (Exception ex) + { + return ProbeResult.Failure("Open.Numeric.Primes.Prime.Numbers.IsPrime", ex); + } +} + +static void OpenNumericPrimesIsPrimeKernel(Index1D index, ArrayView1D data) +{ + ulong value = (ulong)(index + 2); + bool isPrime = Prime.Numbers.IsPrime(value); + data[index] = isPrime ? 1 : 0; +} + static void BigIntegerKernel(Index1D index, ArrayView1D data) { BigInteger value = new BigInteger(data[index]);