From 38df505f4fd6c86d5565e733f3d40e7976c55a8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 4 Oct 2025 22:40:12 +0200 Subject: [PATCH] Implement GPU pow2 heuristics with safe fallbacks --- .../Gpu/PrimeOrderGpuHeuristicsTests.cs | 18 +-- .../Gpu/PrimeOrderGpuHeuristics.cs | 141 +++++++++++++----- PerfectNumbers.Core/PrimeOrderCalculator.cs | 1 - PerfectNumbers.Core/ULongExtensions.cs | 11 ++ todo.txt | 2 - 5 files changed, 123 insertions(+), 50 deletions(-) diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs index 9b0a40c..1bb9335 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs @@ -35,12 +35,12 @@ public class PrimeOrderGpuHeuristicsTests } [Fact] - public void TryPow2Mod_returns_unavailable_for_unmarked_prime() + public void TryPow2Mod_returns_success_for_supported_prime() { GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, 193UL, out ulong remainder); - status.Should().Be(GpuPow2ModStatus.Unavailable); - remainder.Should().Be(0UL); + status.Should().Be(GpuPow2ModStatus.Success); + remainder.Should().Be(2UL); } [Fact] @@ -76,8 +76,8 @@ public class PrimeOrderGpuHeuristicsTests remainder.Should().Be(0UL); GpuPow2ModStatus followUp = PrimeOrderGpuHeuristics.TryPow2Mod(4UL, 193UL, out ulong followUpRemainder); - followUp.Should().Be(GpuPow2ModStatus.Unavailable); - followUpRemainder.Should().Be(0UL); + followUp.Should().Be(GpuPow2ModStatus.Success); + followUpRemainder.Should().Be(16UL); } finally { @@ -108,15 +108,15 @@ public class PrimeOrderGpuHeuristicsTests } [Fact] - public void TryPow2ModBatch_returns_unavailable_and_clears_remainders_for_unmarked_prime() + public void TryPow2ModBatch_returns_success_for_supported_prime() { ulong[] remainders = { 123UL, 456UL }; ulong[] exponents = { 1UL, 2UL }; GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 257UL, remainders); - status.Should().Be(GpuPow2ModStatus.Unavailable); - remainders[0].Should().Be(0UL); - remainders[1].Should().Be(0UL); + status.Should().Be(GpuPow2ModStatus.Success); + remainders[0].Should().Be(2UL); + remainders[1].Should().Be(4UL); } [Fact] diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 8a5643a..433e5aa 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -1,5 +1,9 @@ using System; using System.Collections.Concurrent; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core; +using System.Runtime.InteropServices; using System.Numerics; namespace PerfectNumbers.Core.Gpu; @@ -14,6 +18,7 @@ internal enum GpuPow2ModStatus internal static class PrimeOrderGpuHeuristics { private static readonly ConcurrentDictionary OverflowedPrimes = new(); + private static readonly ConcurrentDictionary, GpuUInt128, ArrayView1D>> Pow2ModKernelCache = new(); private static PrimeOrderGpuCapability s_capability = PrimeOrderGpuCapability.Default; internal static ConcurrentDictionary OverflowRegistry => OverflowedPrimes; @@ -30,34 +35,13 @@ internal static class PrimeOrderGpuHeuristics public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder) { - remainder = 0UL; - if (prime <= 1UL) - { - return GpuPow2ModStatus.Unavailable; - } + Span exponents = stackalloc ulong[1]; + Span remainders = stackalloc ulong[1]; + exponents[0] = exponent; - ConcurrentDictionary overflowRegistry = OverflowedPrimes; - - if (overflowRegistry.ContainsKey(prime)) - { - return GpuPow2ModStatus.Overflow; - } - - PrimeOrderGpuCapability capability = s_capability; - - if (!SupportsPrime(prime, capability)) - { - overflowRegistry[prime] = 0; - return GpuPow2ModStatus.Overflow; - } - - if (!SupportsExponent(exponent, capability)) - { - return GpuPow2ModStatus.Overflow; - } - - // TODO: Implement the GPU heuristic order calculator so this path mirrors PrimeOrderCalculator before enabling device execution. - return GpuPow2ModStatus.Unavailable; + GpuPow2ModStatus status = TryPow2ModBatch(exponents, prime, remainders); + remainder = remainders[0]; + return status; } public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan exponents, ulong prime, Span remainders) @@ -72,7 +56,8 @@ internal static class PrimeOrderGpuHeuristics throw new ArgumentException("Remainder span is shorter than the exponent span.", nameof(remainders)); } - remainders.Slice(0, exponents.Length).Clear(); + Span target = remainders.Slice(0, exponents.Length); + target.Clear(); if (prime <= 1UL) { @@ -102,33 +87,113 @@ internal static class PrimeOrderGpuHeuristics } } - // TODO: Implement the GPU heuristic order calculator so this path mirrors PrimeOrderCalculator before enabling device execution. - return GpuPow2ModStatus.Unavailable; + bool computed = TryComputeOnGpu(exponents, prime, target); + return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable; } private static bool SupportsPrime(ulong prime, PrimeOrderGpuCapability capability) { - return GetBitLength(prime) <= capability.ModulusBits; + return prime.GetBitLength() <= capability.ModulusBits; } private static bool SupportsExponent(ulong exponent, PrimeOrderGpuCapability capability) { - return GetBitLength(exponent) <= capability.ExponentBits; + return exponent.GetBitLength() <= capability.ExponentBits; } - private static int GetBitLength(ulong value) + private static bool TryComputeOnGpu(ReadOnlySpan exponents, ulong prime, Span results) { - if (value == 0UL) + try { - return 0; + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + try + { + using var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + var kernel = GetPow2ModKernel(accelerator); + using var exponentBuffer = accelerator.Allocate1D(exponents.Length); + using var remainderBuffer = accelerator.Allocate1D(exponents.Length); + + exponentBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(exponents), exponents.Length); + remainderBuffer.MemSetToZero(); + + kernel(stream, exponents.Length, exponentBuffer.View, new GpuUInt128(prime), remainderBuffer.View); + + stream.Synchronize(); + + remainderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), results.Length); + return true; + } + finally + { + lease.Dispose(); + } + } + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + { + ComputePow2ModCpu(exponents, prime, results); + return true; + } + } + + private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, Span results) + { + for (int i = 0; i < exponents.Length; i++) + { + results[i] = Pow2ModCpu(exponents[i], prime); + } + } + + private static ulong Pow2ModCpu(ulong exponent, ulong modulus) + { + if (modulus <= 1UL) + { + return 0UL; } - return 64 - BitOperations.LeadingZeroCount(value); + ulong result = 1UL % modulus; + ulong baseValue = 2UL % modulus; + ulong remaining = exponent; + + while (remaining > 0UL) + { + if ((remaining & 1UL) != 0UL) + { + result = (ulong)(((UInt128)result * baseValue) % modulus); + } + + remaining >>= 1; + if (remaining == 0UL) + { + break; + } + + baseValue = (ulong)(((UInt128)baseValue * baseValue) % modulus); + } + + return result; + } + + private static Action, GpuUInt128, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) + { + return Pow2ModKernelCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(Pow2ModKernel); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); + }); + } + + private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, GpuUInt128 modulus, ArrayView1D remainders) + { + ulong exponent = exponents[index]; + GpuUInt128 result = GpuUInt128.Pow2Mod(exponent, modulus); + remainders[index] = result.Low; } internal readonly record struct PrimeOrderGpuCapability(int ModulusBits, int ExponentBits) { public static PrimeOrderGpuCapability Default => new(128, 64); } - - } +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 93d91e0..240d4ed 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -691,7 +691,6 @@ internal static partial class PrimeOrderCalculator private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) { - // TODO: Implement the GPU heuristic path so TryPow2Mod can verify candidates on the device instead of immediately falling back to the CPU Montgomery reducer. GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder); if (status == GpuPow2ModStatus.Success) { diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index f94b2eb..58f2802 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -56,6 +56,17 @@ public static class ULongExtensions return order; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetBitLength(this ulong value) + { + if (value == 0UL) + { + return 0; + } + + return 64 - BitOperations.LeadingZeroCount(value); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MultiplyShiftRight(ulong value, ulong multiplier, int shift) { diff --git a/todo.txt b/todo.txt index ca5be75..19278dc 100644 --- a/todo.txt +++ b/todo.txt @@ -1,3 +1 @@ -- Move the ulong bit-length helper into ULongExtensions and update PrimeOrderGpuHeuristics (and any other callers) to reuse it. -- Implement the GPU heuristics for PrimeOrderCalculator using ILGPU. Reuse the numeric type compatibility probe results to select between EInteger, ERational, GpuUInt128, or any new GPU-specific big-integer representations, guarantee parity with the CPU heuristic algorithm (including 128-bit primes), prefer the shift-first multiply/shift variant to stay within device ranges, and wire the calculator through the accelerator path. - 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.