From 6d82f5a608b178f450dde1d555241a4c7d40d5f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Fri, 17 Oct 2025 04:11:16 +0200 Subject: [PATCH] Enable GPU execution for windowed Montgomery pow2 --- ...MersenneNumberDivisorByDivisorGpuTester.cs | 4 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 8 +- PerfectNumbers.Core/ULongExtensions.cs | 152 ++++++++++++++---- 3 files changed, 130 insertions(+), 34 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 2e77781..909176e 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -546,7 +546,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; + hits[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -765,7 +765,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong exponent = exponents[index]; - results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true); + results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index bd703a7..7ce5c57 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -809,7 +809,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = phi / factor; - if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL) + if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL) { return false; } @@ -1593,7 +1593,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) { - return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL; + return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL; } private static ulong CalculateByDoublingKernel(ulong prime) @@ -1793,7 +1793,7 @@ internal static partial class PrimeOrderGpuHeuristics return 0UL; } - return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false); + return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisorData, keepMontgomery: false); } private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) @@ -1809,7 +1809,7 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { ulong exponent = exponents[index]; - remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false); + remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index fa2ba36..cc71a47 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -1,6 +1,9 @@ +using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; +using ILGPU; using ILGPU.Algorithms; +using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -409,13 +412,32 @@ public static class ULongExtensions return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); } + if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult)) + { + return gpuResult; + } + + return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong Pow2MontgomeryModWindowedKernel(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + ulong modulus = divisor.Modulus; + if (exponent == 0UL) + { + return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + } + + if (exponent <= Pow2WindowFallbackThreshold) + { + return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + } + int bitLength = GetPortableBitLength(exponent); int windowSize = GetWindowSize(bitLength); - int oddPowerCount = 1 << (windowSize - 1); - ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - ulong[] oddPowers = InitializeMontgomeryOddPowersGpu(divisor, modulus, nPrime, oddPowerCount); int index = bitLength - 1; while (index >= 0) @@ -446,19 +468,113 @@ public static class ULongExtensions ulong mask = (1UL << windowLength) - 1UL; ulong windowValue = (exponent >> windowStart) & mask; - int tableIndex = (int)((windowValue - 1UL) >> 1); - ulong multiplier = oddPowers[tableIndex]; + ulong multiplier = ComputeMontgomeryOddPower(windowValue, divisor, modulus, nPrime); result = result.MontgomeryMultiply(multiplier, modulus, nPrime); index = windowStart - 1; } - if (keepMontgomery) + if (!keepMontgomery) { - return result; + result = result.MontgomeryMultiply(1UL, modulus, nPrime); } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeMontgomeryOddPower(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime) + { + ulong baseValue = divisor.MontgomeryTwo; + ulong power = divisor.MontgomeryOne; + ulong remaining = exponent; + + while (remaining != 0UL) + { + if ((remaining & 1UL) != 0UL) + { + power = power.MontgomeryMultiply(baseValue, modulus, nPrime); + } + + remaining >>= 1; + if (remaining == 0UL) + { + break; + } + + baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); + } + + return power; + } + + private static class Pow2MontgomeryGpuExecutor + { + private static readonly ConcurrentDictionary, MontgomeryDivisorData, byte, ArrayView1D>> KernelCache = new(); + + public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result) + { + result = 0UL; + + GpuKernelLease lease; + try + { + lease = GpuKernelPool.GetKernel(useGpuOrder: true); + } + catch (Exception) + { + return false; + } + + var execution = lease.EnterExecutionScope(); + try + { + Accelerator accelerator = lease.Accelerator; + if (accelerator.AcceleratorType == AcceleratorType.CPU) + { + return false; + } + + var kernel = KernelCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernel); + var launcher = KernelUtil.GetKernel(loaded); + return launcher.CreateLauncherDelegate, MontgomeryDivisorData, byte, ArrayView1D>>(); + }); + + using var exponentBuffer = accelerator.Allocate1D(1); + using var resultBuffer = accelerator.Allocate1D(1); + + exponentBuffer.View.CopyFromCPU(ref exponent, 1); + resultBuffer.MemSetToZero(); + + byte keepFlag = keepMontgomery ? (byte)1 : (byte)0; + AcceleratorStream stream = lease.Stream; + kernel(stream, 1, exponentBuffer.View, divisor, keepFlag, resultBuffer.View); + stream.Synchronize(); + + resultBuffer.View.CopyToCPU(ref result, 1); + return true; + } + catch (AcceleratorException) + { + return false; + } + catch (NotSupportedException) + { + return false; + } + finally + { + execution.Dispose(); + lease.Dispose(); + } + } + + private static void Pow2MontgomeryKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, byte keepMontgomery, ArrayView1D results) + { + results[index] = Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0); + } } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -549,26 +665,6 @@ public static class ULongExtensions } } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) - { - ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount]; - oddPowers[0] = divisor.MontgomeryTwo; - if (oddPowerCount == 1) - { - return oddPowers; - } - - ulong square = divisor.MontgomeryTwoSquared; - for (int i = 1; i < oddPowerCount; i++) - { - ulong previous = oddPowers[i - 1]; - oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); - } - - return oddPowers; - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) {