diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 6b022aa..cf03ebb 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -37,9 +37,16 @@ public sealed class KernelContainer ResidueAutomatonArgs, ArrayView, ArrayView1D>? IncrementalOrder; public Action, ArrayView1D>? Pow2ModOrder; + public Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>? SmallPrimeFactor; // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length. public MemoryBuffer1D? SmallCycles; + public MemoryBuffer1D? SmallPrimeFactorsPrimes; + public MemoryBuffer1D? SmallPrimeFactorsSquares; + public MemoryBuffer1D? SmallPrimeFactorPrimeSlots; + public MemoryBuffer1D? SmallPrimeFactorExponentSlots; + public MemoryBuffer1D? SmallPrimeFactorCountSlot; + public MemoryBuffer1D? SmallPrimeFactorRemainingSlot; public MemoryBuffer1D? SmallPrimesLastOne; public MemoryBuffer1D? SmallPrimesPow2LastOne; public MemoryBuffer1D? SmallPrimesLastSeven; @@ -75,6 +82,54 @@ public readonly struct ResiduePrimeViews( public readonly ArrayView1D LastSevenPow2 = lastSevenPow2; } +public readonly struct SmallPrimeFactorTables +{ + public readonly MemoryBuffer1D Primes; + public readonly MemoryBuffer1D Squares; + public readonly int Count; + + public SmallPrimeFactorTables( + MemoryBuffer1D primes, + MemoryBuffer1D squares, + int count) + { + Primes = primes; + Squares = squares; + Count = count; + } + + public ArrayView1D PrimesView => Primes.View; + public ArrayView1D SquaresView => Squares.View; +} + +public readonly struct SmallPrimeFactorScratch +{ + public readonly MemoryBuffer1D PrimeSlots; + public readonly MemoryBuffer1D ExponentSlots; + public readonly MemoryBuffer1D CountSlot; + public readonly MemoryBuffer1D RemainingSlot; + + public SmallPrimeFactorScratch( + MemoryBuffer1D primeSlots, + MemoryBuffer1D exponentSlots, + MemoryBuffer1D countSlot, + MemoryBuffer1D remainingSlot) + { + PrimeSlots = primeSlots; + ExponentSlots = exponentSlots; + CountSlot = countSlot; + RemainingSlot = remainingSlot; + } + + public void Clear() + { + PrimeSlots.MemSetToZero(); + ExponentSlots.MemSetToZero(); + CountSlot.MemSetToZero(); + RemainingSlot.MemSetToZero(); + } +} + public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels) { private bool disposedValue; @@ -189,6 +244,82 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } } + public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernel + { + get + { + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.SmallPrimeFactor, () => + { + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + }); + } + } + } + + private static void SmallPrimeFactorKernelScan( + Index1D index, + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D primeSlots, + ArrayView1D exponentSlots, + ArrayView1D countSlot, + ArrayView1D remainingSlot) + { + if (index != 0) + { + return; + } + + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + ulong remaining = value; + int count = 0; + int capacity = (int)primeSlots.Length; + + for (int i = 0; i < primeCount && remaining > 1UL && count < capacity; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + primeSlots[count] = primeValue; + exponentSlots[count] = exponent; + count++; + } + + countSlot[0] = count; + remainingSlot[0] = remaining; + } + // TODO: Plumb the small-cycles device buffer into all kernels that can benefit // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. @@ -779,6 +910,75 @@ public class GpuKernelPool } } + public static SmallPrimeFactorTables EnsureSmallPrimeFactorTables(Accelerator accelerator) + { + var kernels = GetKernels(accelerator); + if (kernels.SmallPrimeFactorsPrimes is { } primeBuffer && kernels.SmallPrimeFactorsSquares is { } squareBuffer) + { + return new SmallPrimeFactorTables(primeBuffer, squareBuffer, (int)primeBuffer.Length); + } + + lock (kernels) + { + if (kernels.SmallPrimeFactorsPrimes is { } existingPrimes && kernels.SmallPrimeFactorsSquares is { } existingSquares) + { + return new SmallPrimeFactorTables(existingPrimes, existingSquares, (int)existingPrimes.Length); + } + + var hostPrimes = PrimesGenerator.SmallPrimes; + var hostSquares = PrimesGenerator.SmallPrimesPow2; + + var devicePrimes = accelerator.Allocate1D(hostPrimes.Length); + devicePrimes.View.CopyFromCPU(hostPrimes); + var deviceSquares = accelerator.Allocate1D(hostSquares.Length); + deviceSquares.View.CopyFromCPU(hostSquares); + + kernels.SmallPrimeFactorsPrimes = devicePrimes; + kernels.SmallPrimeFactorsSquares = deviceSquares; + + return new SmallPrimeFactorTables(devicePrimes, deviceSquares, hostPrimes.Length); + } + } + + public static SmallPrimeFactorScratch EnsureSmallPrimeFactorScratch(Accelerator accelerator, int slotCount) + { + var kernels = GetKernels(accelerator); + lock (kernels) + { + MemoryBuffer1D? primeSlots = kernels.SmallPrimeFactorPrimeSlots; + if (primeSlots is null || primeSlots.Length < slotCount) + { + primeSlots?.Dispose(); + primeSlots = accelerator.Allocate1D(slotCount); + kernels.SmallPrimeFactorPrimeSlots = primeSlots; + } + + MemoryBuffer1D? exponentSlots = kernels.SmallPrimeFactorExponentSlots; + if (exponentSlots is null || exponentSlots.Length < slotCount) + { + exponentSlots?.Dispose(); + exponentSlots = accelerator.Allocate1D(slotCount); + kernels.SmallPrimeFactorExponentSlots = exponentSlots; + } + + MemoryBuffer1D? countSlot = kernels.SmallPrimeFactorCountSlot; + if (countSlot is null) + { + countSlot = accelerator.Allocate1D(1); + kernels.SmallPrimeFactorCountSlot = countSlot; + } + + MemoryBuffer1D? remainingSlot = kernels.SmallPrimeFactorRemainingSlot; + if (remainingSlot is null) + { + remainingSlot = accelerator.Allocate1D(1); + kernels.SmallPrimeFactorRemainingSlot = remainingSlot; + } + + return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!); + } + } + public static GpuKernelLease GetKernel(bool useGpuOrder) { var limiter = GpuPrimeWorkLimiter.Acquire(); 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/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index cedaff1..f21d726 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,3 +1,8 @@ +using System; +using System.Runtime.InteropServices; +using ILGPU; +using ILGPU.Runtime; +using ILGPU.Runtime.OpenCL; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -16,30 +21,88 @@ internal static partial class PrimeOrderCalculator Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); primeBuffer.Clear(); exponentBuffer.Clear(); + remaining = value; - if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) + try + { + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + try + { + var execution = lease.EnterExecutionScope(); + try + { + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator); + SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); + scratch.Clear(); + + var kernel = lease.SmallPrimeFactorKernel; + kernel( + stream, + 1, + value, + limit, + tables.PrimesView, + tables.SquaresView, + tables.Count, + scratch.PrimeSlots.View, + scratch.ExponentSlots.View, + scratch.CountSlot.View, + scratch.RemainingSlot.View); + + stream.Synchronize(); + + int factorCount = 0; + scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); + factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); + + if (factorCount > 0) + { + scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount); + scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount); + } + scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); + + 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; + } + finally + { + execution.Dispose(); + } + } + finally + { + lease.Dispose(); + } + } + catch (CLException) { remaining = value; - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); return false; } - - remaining = gpuRemaining; - for (int i = 0; i < factorCount; i++) + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - if (primeValue == 0UL || exponent == 0) - { - continue; - } - - counts[primeValue] = exponent; + remaining = value; + return false; + } + finally + { + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); } - - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - return true; } } 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) {