From 6458f6740d3173591334f0852389e8ee120ca6b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 20:30:44 +0100 Subject: [PATCH 1/7] Add GPU prime tester overload for preallocated buffers --- .../Gpu/PrimeTesterGpuTests.cs | 30 +++++ PerfectNumbers.Core/PrimeTester.cs | 114 +++++++++++++++--- 2 files changed, 127 insertions(+), 17 deletions(-) diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs index 5e8faf9..0cce423 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs @@ -1,4 +1,6 @@ using FluentAssertions; +using ILGPU; +using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; using Xunit; @@ -27,6 +29,34 @@ public class PrimeTesterGpuTests } } + [Fact] + public void IsPrimeBatchGpu_with_preallocated_resources_marks_expected_values() + { + GpuContextPool.ForceCpu = false; + using var context = Context.CreateDefault(); + using var accelerator = context.GetPreferredDevice(false).CreateAccelerator(context); + + try + { + ulong[] values = [31UL, 33UL, 37UL, 39UL]; + byte[] results = new byte[values.Length]; + + using var inputBuffer = accelerator.Allocate1D(values.Length); + using var outputBuffer = accelerator.Allocate1D(values.Length); + ulong[] staging = new ulong[values.Length]; + + PrimeTester.IsPrimeBatchGpu(values, results, accelerator, inputBuffer, outputBuffer, staging); + + results.Should().Equal(new byte[] { 1, 0, 1, 0 }); + } + finally + { + PrimeTester.ClearGpuCaches(accelerator); + GpuContextPool.ForceCpu = false; + GpuContextPool.DisposeAll(); + } + } + [Fact] public void IsPrimeGpu_rejects_composites() { diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index dd15c4e..ba46038 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -3,6 +3,7 @@ using System.Buffers; using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; using System.Threading; using ILGPU; using ILGPU.Algorithms; @@ -112,7 +113,6 @@ public sealed class PrimeTester var gpu = PrimeTesterGpuContextPool.Rent(); var accelerator = gpu.Accelerator; var state = GpuKernelState.GetOrCreate(accelerator); - int totalLength = values.Length; int batchSize = Math.Max(1, GpuBatchSize); lock (gpu.ExecutionLock) @@ -123,30 +123,110 @@ public sealed class PrimeTester ArrayPool pool = ThreadStaticPools.UlongPool; ulong[] temp = pool.Rent(batchSize); - int pos = 0; - while (pos < totalLength) + try { - int remaining = totalLength - pos; - int count = remaining > batchSize ? batchSize : remaining; - - values.Slice(pos, count).CopyTo(temp); - input.View.CopyFromCPU(ref temp[0], count); - - state.Kernel(count, input.View, state.DevicePrimes.View, output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref results[pos], count); - - pos += count; + IsPrimeBatchGpu(values, results, accelerator, input, output, temp); + } + finally + { + pool.Return(temp, clearArray: false); + state.ReturnScratch(scratch); } - - pool.Return(temp, clearArray: false); - state.ReturnScratch(scratch); } gpu.Dispose(); limiter.Dispose(); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void IsPrimeBatchGpu( + ReadOnlySpan values, + Span results, + Accelerator accelerator, + MemoryBuffer1D deviceInput, + MemoryBuffer1D deviceOutput, + Span stagingBuffer) + { + if (accelerator is null) + { + throw new ArgumentNullException(nameof(accelerator)); + } + + if (deviceInput is null) + { + throw new ArgumentNullException(nameof(deviceInput)); + } + + if (deviceOutput is null) + { + throw new ArgumentNullException(nameof(deviceOutput)); + } + + if (results.Length < values.Length) + { + throw new ArgumentException("The results span must be at least as long as the input values.", nameof(results)); + } + + if (values.IsEmpty) + { + return; + } + + if (stagingBuffer.IsEmpty) + { + throw new ArgumentException("The staging buffer must contain at least one element.", nameof(stagingBuffer)); + } + + long inputLength = deviceInput.Length; + long outputLength = deviceOutput.Length; + if (inputLength <= 0) + { + throw new ArgumentException("The device input buffer must contain at least one element.", nameof(deviceInput)); + } + + if (outputLength <= 0) + { + throw new ArgumentException("The device output buffer must contain at least one element.", nameof(deviceOutput)); + } + + long stagingLength = stagingBuffer.Length; + long capacityLong = Math.Min(Math.Min(inputLength, outputLength), stagingLength); + if (capacityLong <= 0) + { + throw new ArgumentException("The provided buffers must be able to hold at least one element."); + } + + int capacity = (int)Math.Min(capacityLong, int.MaxValue); + + var state = GpuKernelState.GetOrCreate(accelerator); + var inputView = deviceInput.View; + var outputView = deviceOutput.View; + var primesView = state.DevicePrimes.View; + var stagingSpan = stagingBuffer; + ref ulong stagingRef = ref MemoryMarshal.GetReference(stagingSpan); + + int totalLength = values.Length; + int pos = 0; + while (pos < totalLength) + { + int remaining = totalLength - pos; + int count = remaining > capacity ? capacity : remaining; + + var stagingSlice = stagingSpan.Slice(0, count); + values.Slice(pos, count).CopyTo(stagingSlice); + + var inputSlice = inputView.SubView(0, count); + inputSlice.CopyFromCPU(ref stagingRef, count); + + var outputSlice = outputView.SubView(0, count); + state.Kernel(count, inputSlice, primesView, outputSlice); + accelerator.Synchronize(); + outputSlice.CopyToCPU(ref results[pos], count); + + pos += count; + } + } + internal static class PrimeTesterGpuContextPool { internal sealed class PooledContext : IDisposable From ba439f0f0fba496a02e0133645710d9acc76eacb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 20:58:21 +0100 Subject: [PATCH 2/7] Reuse GPU context for nested partial factoring --- PerfectNumbers.Core/Gpu/GpuKernelPool.cs | 16 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 198 ++++++++++-------- .../PrimeOrderCalculator.Cpu.cs | 134 ++++++++++-- .../PrimeOrderCalculator.Gpu.cs | 68 +++--- 4 files changed, 270 insertions(+), 146 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index c16d9ca..2060b4f 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -46,6 +46,7 @@ public sealed class KernelContainer public MemoryBuffer1D? SmallPrimeFactorExponentSlots; public MemoryBuffer1D? SmallPrimeFactorCountSlot; public MemoryBuffer1D? SmallPrimeFactorRemainingSlot; + public MemoryBuffer1D? SmallPrimeFactorFullyFactoredSlot; public MemoryBuffer1D? SmallPrimesLastOne; public MemoryBuffer1D? SmallPrimesPow2LastOne; public MemoryBuffer1D? SmallPrimesLastSeven; @@ -107,17 +108,20 @@ public readonly struct SmallPrimeFactorScratch public readonly MemoryBuffer1D ExponentSlots; public readonly MemoryBuffer1D CountSlot; public readonly MemoryBuffer1D RemainingSlot; + public readonly MemoryBuffer1D FullyFactoredSlot; public SmallPrimeFactorScratch( MemoryBuffer1D primeSlots, MemoryBuffer1D exponentSlots, MemoryBuffer1D countSlot, - MemoryBuffer1D remainingSlot) + MemoryBuffer1D remainingSlot, + MemoryBuffer1D fullyFactoredSlot) { PrimeSlots = primeSlots; ExponentSlots = exponentSlots; CountSlot = countSlot; RemainingSlot = remainingSlot; + FullyFactoredSlot = fullyFactoredSlot; } public void Clear() @@ -126,6 +130,7 @@ public readonly struct SmallPrimeFactorScratch ExponentSlots.MemSetToZero(); CountSlot.MemSetToZero(); RemainingSlot.MemSetToZero(); + FullyFactoredSlot.MemSetToZero(); } } @@ -273,7 +278,14 @@ public class GpuKernelPool kernels.SmallPrimeFactorRemainingSlot = remainingSlot; } - return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!); + MemoryBuffer1D? fullyFactoredSlot = kernels.SmallPrimeFactorFullyFactoredSlot; + if (fullyFactoredSlot is null) + { + fullyFactoredSlot = accelerator.Allocate1D(1); + kernels.SmallPrimeFactorFullyFactoredSlot = fullyFactoredSlot; + } + + return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!, fullyFactoredSlot!); } } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 7b17b48..04ed0ce 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -151,34 +151,52 @@ internal static partial class PrimeOrderGpuHeuristics primeTargets.Clear(); exponentTargets.Clear(); - if (!TryLaunchPartialFactorKernel( - value, - limit, - primeTargets, - exponentTargets, - out int extracted, - out ulong leftover, - out bool kernelFullyFactored)) + try { - return false; + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + try + { + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, primeTargets.Length); + return TryPartialFactor( + accelerator, + stream, + scratch, + value, + limit, + primeTargets, + exponentTargets, + out factorCount, + out remaining, + out fullyFactored); + } + finally + { + execution.Dispose(); + lease.Dispose(); + } + } + catch (CLException ex) + { + Console.WriteLine($"GPU ERROR ({ex.Error}): {ex.Message}"); + } + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + { + Console.WriteLine($"GPU ERROR: {ex.Message}"); } - 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; + factorCount = 0; + remaining = value; + fullyFactored = false; + return false; } - private static bool TryLaunchPartialFactorKernel( + public static bool TryPartialFactor( + Accelerator accelerator, + AcceleratorStream stream, + SmallPrimeFactorScratch scratch, ulong value, uint limit, Span primeTargets, @@ -191,83 +209,99 @@ internal static partial class PrimeOrderGpuHeuristics remaining = value; fullyFactored = false; + primeTargets.Clear(); + exponentTargets.Clear(); + try { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - var execution = lease.EnterExecutionScope(); - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - - var kernel = GetPartialFactorKernel(accelerator); - SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); - - // TODO: We should create / reallocate these buffer only once or if the new required length is bigger than capacity. - var factorBuffer = accelerator.Allocate1D(primeTargets.Length); - var exponentBuffer = accelerator.Allocate1D(exponentTargets.Length); - var countBuffer = accelerator.Allocate1D(1); - var remainingBuffer = accelerator.Allocate1D(1); - var fullyFactoredBuffer = accelerator.Allocate1D(1); - - // TODO: There is no need to clear these buffers because the kernel will always assign values within the required bounds. - // factorBuffer.MemSetToZero(); - // exponentBuffer.MemSetToZero(); - // countBuffer.MemSetToZero(); - // remainingBuffer.MemSetToZero(); - // fullyFactoredBuffer.MemSetToZero(); - - kernel( + return TryPartialFactorCore( + accelerator, stream, - 1, - cache.Primes!.View, - cache.Squares!.View, - cache.Count, - primeTargets.Length, + scratch, 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; - - // TODO: These buffers shouldn't be disposed but rather reused accross call with the assigned accelerator. - factorBuffer.Dispose(); - exponentBuffer.Dispose(); - countBuffer.Dispose(); - remainingBuffer.Dispose(); - fullyFactoredBuffer.Dispose(); - execution.Dispose(); - lease.Dispose(); - return true; + primeTargets, + exponentTargets, + out factorCount, + out remaining, + out fullyFactored); } catch (CLException ex) { Console.WriteLine($"GPU ERROR ({ex.Error}): {ex.Message}"); - factorCount = 0; - remaining = value; - fullyFactored = false; - return false; } catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) { Console.WriteLine($"GPU ERROR: {ex.Message}"); + } + + factorCount = 0; + remaining = value; + fullyFactored = false; + return false; + } + + private static bool TryPartialFactorCore( + Accelerator accelerator, + AcceleratorStream stream, + SmallPrimeFactorScratch scratch, + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining, + out bool fullyFactored) + { + factorCount = 0; + remaining = value; + fullyFactored = false; + + scratch.Clear(); + + var kernel = GetPartialFactorKernel(accelerator); + SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); + + kernel( + stream, + 1, + cache.Primes!.View, + cache.Squares!.View, + cache.Count, + primeTargets.Length, + value, + limit, + scratch.PrimeSlots.View, + scratch.ExponentSlots.View, + scratch.CountSlot.View, + scratch.RemainingSlot.View, + scratch.FullyFactoredSlot.View); + + stream.Synchronize(); + + scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + if (factorCount > capacity) + { factorCount = 0; remaining = value; fullyFactored = false; return false; } + + if (factorCount > 0) + { + scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeTargets), factorCount); + scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentTargets), factorCount); + } + + scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); + + byte fullyFactoredFlag = 0; + scratch.FullyFactoredSlot.View.CopyToCPU(ref fullyFactoredFlag, 1); + bool kernelFullyFactored = fullyFactoredFlag != 0; + fullyFactored = kernelFullyFactored && remaining == 1UL; + return true; } public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData divisorData) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 9b83a9c..70e8b6c 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -3,9 +3,10 @@ using System.Buffers; using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; -using System.Runtime.InteropServices; using ILGPU; +using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; +using System.Runtime.InteropServices; namespace PerfectNumbers.Core; @@ -43,6 +44,43 @@ internal static partial class PrimeOrderCalculator } } + private sealed class PartialFactorGpuContext : IDisposable + { + private readonly GpuKernelLease _lease; + private readonly GpuKernelLease.ExecutionScope _executionScope; + + private PartialFactorGpuContext(GpuKernelLease lease, GpuKernelLease.ExecutionScope executionScope, Accelerator accelerator, AcceleratorStream stream, SmallPrimeFactorScratch scratch) + { + _lease = lease; + _executionScope = executionScope; + Accelerator = accelerator; + Stream = stream; + Scratch = scratch; + } + + public Accelerator Accelerator { get; } + + public AcceleratorStream Stream { get; } + + public SmallPrimeFactorScratch Scratch { get; } + + public static PartialFactorGpuContext Rent() + { + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); + return new PartialFactorGpuContext(lease, execution, accelerator, stream, scratch); + } + + public void Dispose() + { + _executionScope.Dispose(); + _lease.Dispose(); + } + } + private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. @@ -58,19 +96,38 @@ internal static partial class PrimeOrderCalculator return gpuOrder; } - PartialFactorResult phiFactors = PartialFactor(phi, config); + PartialFactorGpuContext? gpuContext = null; + bool disposeContext = false; - ulong result; - if (phiFactors.Factors is null) + if (IsGpuHeuristicDevice) { - result = CalculateByFactorizationCpu(prime, divisorData); + gpuContext = PartialFactorGpuContext.Rent(); + disposeContext = true; + } + + try + { + PartialFactorResult phiFactors = PartialFactor(phi, config, gpuContext); + + ulong result; + if (phiFactors.Factors is null) + { + result = CalculateByFactorizationCpu(prime, divisorData); + phiFactors.Dispose(); + return result; + } + + result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors, gpuContext); phiFactors.Dispose(); return result; } - - result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); - phiFactors.Dispose(); - return result; + finally + { + if (disposeContext && gpuContext is not null) + { + gpuContext.Dispose(); + } + } } private static ulong RunHeuristicPipelineCpu( @@ -79,7 +136,8 @@ internal static partial class PrimeOrderCalculator in PrimeOrderSearchConfig config, in MontgomeryDivisorData divisorData, ulong phi, - PartialFactorResult phiFactors) + PartialFactorResult phiFactors, + PartialFactorGpuContext? gpuContext) { if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) { @@ -89,7 +147,7 @@ internal static partial class PrimeOrderCalculator ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, gpuContext, out PartialFactorResult? orderFactors)) { return candidateOrder; } @@ -100,7 +158,7 @@ internal static partial class PrimeOrderCalculator return CalculateByFactorizationCpu(prime, divisorData); } - if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) + if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, gpuContext, out ulong order)) { return order; } @@ -330,6 +388,7 @@ internal static partial class PrimeOrderCalculator ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, + PartialFactorGpuContext? gpuContext, out PartialFactorResult? reusableFactorization) { reusableFactorization = null; @@ -344,7 +403,7 @@ internal static partial class PrimeOrderCalculator return false; } - PartialFactorResult factorization = PartialFactor(order, config); + PartialFactorResult factorization = PartialFactor(order, config, gpuContext); if (factorization.Factors is null) { factorization.Dispose(); @@ -492,6 +551,7 @@ internal static partial class PrimeOrderCalculator in PrimeOrderSearchConfig config, PartialFactorResult phiFactors, PartialFactorResult? cachedOrderFactors, + PartialFactorGpuContext? gpuContext, out ulong result) { result = 0UL; @@ -502,7 +562,7 @@ internal static partial class PrimeOrderCalculator } // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactor(order, config); + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactor(order, config, gpuContext); try { if (orderFactors.Factors is null) @@ -618,7 +678,7 @@ internal static partial class PrimeOrderCalculator bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); powStepperInitialized = true; - if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, gpuContext, ref powUsed, powBudget)) { if (gpuPool is not null) { @@ -655,7 +715,7 @@ internal static partial class PrimeOrderCalculator continue; } - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, gpuContext, ref powUsed, powBudget)) { continue; } @@ -858,9 +918,9 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) + private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, PartialFactorGpuContext? gpuContext, ref int powUsed, int powBudget) { - PartialFactorResult factorization = PartialFactor(candidate, config); + PartialFactorResult factorization = PartialFactor(candidate, config, gpuContext); try { if (factorization.Factors is null) @@ -1024,7 +1084,7 @@ internal static partial class PrimeOrderCalculator // private static ulong _partialFactorPendingHits; // private static ulong _partialFactorCofactorHits; - private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config, PartialFactorGpuContext? gpuContext = null) { if (value <= 1UL) { @@ -1033,6 +1093,17 @@ internal static partial class PrimeOrderCalculator const int FactorSlotCount = GpuSmallPrimeFactorSlots; + PartialFactorGpuContext? localContext = gpuContext; + bool disposeContext = false; + + if (IsGpuHeuristicDevice && localContext is null) + { + localContext = PartialFactorGpuContext.Rent(); + disposeContext = true; + } + + try + { // stackalloc is faster than pooling // ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); // int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); @@ -1052,7 +1123,21 @@ internal static partial class PrimeOrderCalculator ulong remaining = value; bool gpuFactored = false; - if (IsGpuHeuristicDevice) + if (IsGpuHeuristicDevice && localContext is not null) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor( + localContext.Accelerator, + localContext.Stream, + localContext.Scratch, + value, + limit, + primeSlots, + exponentSlots, + out factorCount, + out remaining, + out _); + } + else if (IsGpuHeuristicDevice) { gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); } @@ -1329,7 +1414,14 @@ internal static partial class PrimeOrderCalculator ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); return result; } - + finally + { + if (disposeContext && localContext is not null) + { + localContext.Dispose(); + } + } + } private static bool TryPopulateSmallPrimeFactorsCpu( ulong value, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 4cb801c..8e6caf8 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,5 +1,4 @@ using System; -using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; using ILGPU.Runtime.OpenCL; @@ -32,52 +31,39 @@ internal static partial class PrimeOrderCalculator 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); + bool success = PrimeOrderGpuHeuristics.TryPartialFactor( + accelerator, + stream, + scratch, + value, + limit, + primeBuffer, + exponentBuffer, + out factorCount, + out remaining, + out _); - stream.Synchronize(); + if (!success) + { + return false; + } - factorCount = 0; - scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); - factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots); + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + // This will never happen in production code + // if (primeValue == 0UL || exponent == 0) + // { + // continue; + // } - if (factorCount > 0) - { - scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount); - scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount); - } + counts.Add(primeValue, exponent); + } - scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); - - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - // This will never happen in production code - // if (primeValue == 0UL || exponent == 0) - // { - // continue; - // } - - counts.Add(primeValue, exponent); - } - - return true; + return true; } finally { From 794e092e57ae1a3ec222c204b5c19afb7499970b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 21:46:11 +0100 Subject: [PATCH 3/7] Refactor GPU partial factoring context handling --- .../Gpu/PrimeOrderGpuHeuristics.cs | 100 ++--------- .../PrimeOrderCalculator.Cpu.cs | 164 +++++++++--------- .../PrimeOrderCalculator.Gpu.cs | 75 ++++---- PerfectNumbers.Core/PrimeTester.cs | 46 +---- 4 files changed, 129 insertions(+), 256 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 04ed0ce..45ef028 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -130,68 +130,6 @@ internal static partial class PrimeOrderGpuHeuristics }); } - 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(); - - try - { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - var execution = lease.EnterExecutionScope(); - try - { - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, primeTargets.Length); - return TryPartialFactor( - accelerator, - stream, - scratch, - value, - limit, - primeTargets, - exponentTargets, - out factorCount, - out remaining, - out fullyFactored); - } - finally - { - execution.Dispose(); - lease.Dispose(); - } - } - catch (CLException ex) - { - Console.WriteLine($"GPU ERROR ({ex.Error}): {ex.Message}"); - } - catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) - { - Console.WriteLine($"GPU ERROR: {ex.Message}"); - } - - factorCount = 0; - remaining = value; - fullyFactored = false; - return false; - } public static bool TryPartialFactor( Accelerator accelerator, @@ -212,33 +150,17 @@ internal static partial class PrimeOrderGpuHeuristics primeTargets.Clear(); exponentTargets.Clear(); - try - { - return TryPartialFactorCore( - accelerator, - stream, - scratch, - value, - limit, - primeTargets, - exponentTargets, - out factorCount, - out remaining, - out fullyFactored); - } - catch (CLException ex) - { - Console.WriteLine($"GPU ERROR ({ex.Error}): {ex.Message}"); - } - catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) - { - Console.WriteLine($"GPU ERROR: {ex.Message}"); - } - - factorCount = 0; - remaining = value; - fullyFactored = false; - return false; + return TryPartialFactorCore( + accelerator, + stream, + scratch, + value, + limit, + primeTargets, + exponentTargets, + out factorCount, + out remaining, + out fullyFactored); } private static bool TryPartialFactorCore( diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 70e8b6c..4e21cf3 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -48,34 +48,79 @@ internal static partial class PrimeOrderCalculator { private readonly GpuKernelLease _lease; private readonly GpuKernelLease.ExecutionScope _executionScope; + private readonly PrimeTester.KernelState _primeState; + private readonly PrimeTester.KernelState.ScratchBuffers _primeScratch; + private readonly ulong[] _primeHostBuffer; + private readonly byte[] _primeFlagBuffer; + private readonly Action, ArrayView, ArrayView> _primeKernel; + private readonly ArrayView1D _primeDevicePrimes; - private PartialFactorGpuContext(GpuKernelLease lease, GpuKernelLease.ExecutionScope executionScope, Accelerator accelerator, AcceleratorStream stream, SmallPrimeFactorScratch scratch) + private PartialFactorGpuContext( + GpuKernelLease lease, + GpuKernelLease.ExecutionScope executionScope, + SmallPrimeFactorScratch scratch, + PrimeTester.KernelState primeState, + PrimeTester.KernelState.ScratchBuffers primeScratch, + ulong[] primeHostBuffer, + byte[] primeFlagBuffer) { _lease = lease; _executionScope = executionScope; - Accelerator = accelerator; - Stream = stream; Scratch = scratch; + _primeState = primeState; + _primeScratch = primeScratch; + _primeHostBuffer = primeHostBuffer; + _primeFlagBuffer = primeFlagBuffer; + _primeKernel = primeState.Kernel; + _primeDevicePrimes = primeState.DevicePrimes.View; } - public Accelerator Accelerator { get; } + public Accelerator Accelerator => _lease.Accelerator; - public AcceleratorStream Stream { get; } + public AcceleratorStream Stream => _lease.Stream; public SmallPrimeFactorScratch Scratch { get; } + public Span PrimeHostSpan => _primeHostBuffer.AsSpan(0, _primeScratch.Capacity); + + public Span PrimeFlagSpan => _primeFlagBuffer.AsSpan(0, _primeScratch.Capacity); + + public MemoryBuffer1D PrimeDeviceInput => _primeScratch.Input; + + public MemoryBuffer1D PrimeDeviceOutput => _primeScratch.Output; + public static PartialFactorGpuContext Rent() { var lease = GpuKernelPool.GetKernel(useGpuOrder: true); var execution = lease.EnterExecutionScope(); Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); - return new PartialFactorGpuContext(lease, execution, accelerator, stream, scratch); + PrimeTester.KernelState primeState = PrimeTester.GpuKernelState.GetOrCreate(accelerator); + PrimeTester.KernelState.ScratchBuffers primeScratch = primeState.RentScratch(1, accelerator); + ulong[] primeHostBuffer = ThreadStaticPools.UlongPool.Rent(primeScratch.Capacity); + byte[] primeFlagBuffer = ThreadStaticPools.BytePool.Rent(primeScratch.Capacity); + return new PartialFactorGpuContext(lease, execution, scratch, primeState, primeScratch, primeHostBuffer, primeFlagBuffer); + } + + public bool CheckPrime(ulong value) + { + Span host = PrimeHostSpan; + host[0] = value; + ArrayView1D inputSlice = PrimeDeviceInput.View.SubView(0, 1); + inputSlice.CopyFromCPU(ref host[0], 1); + ArrayView1D outputSlice = PrimeDeviceOutput.View.SubView(0, 1); + _primeKernel(1, inputSlice.AsContiguous(), _primeDevicePrimes.AsContiguous(), outputSlice.AsContiguous()); + Accelerator.Synchronize(); + Span flags = PrimeFlagSpan; + outputSlice.CopyToCPU(ref flags[0], 1); + return flags[0] != 0; } public void Dispose() { + _primeState.ReturnScratch(_primeScratch); + ThreadStaticPools.UlongPool.Return(_primeHostBuffer, clearArray: false); + ThreadStaticPools.BytePool.Return(_primeFlagBuffer, clearArray: false); _executionScope.Dispose(); _lease.Dispose(); } @@ -96,38 +141,21 @@ internal static partial class PrimeOrderCalculator return gpuOrder; } - PartialFactorGpuContext? gpuContext = null; - bool disposeContext = false; + PartialFactorGpuContext gpuContext = PartialFactorGpuContext.Rent(); + PartialFactorResult phiFactors = PartialFactor(phi, config, gpuContext); - if (IsGpuHeuristicDevice) + if (phiFactors.Factors is null) { - gpuContext = PartialFactorGpuContext.Rent(); - disposeContext = true; - } - - try - { - PartialFactorResult phiFactors = PartialFactor(phi, config, gpuContext); - - ulong result; - if (phiFactors.Factors is null) - { - result = CalculateByFactorizationCpu(prime, divisorData); - phiFactors.Dispose(); - return result; - } - - result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors, gpuContext); + ulong fallback = CalculateByFactorizationCpu(prime, divisorData); phiFactors.Dispose(); - return result; - } - finally - { - if (disposeContext && gpuContext is not null) - { - gpuContext.Dispose(); - } + gpuContext.Dispose(); + return fallback; } + + ulong result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors, gpuContext); + phiFactors.Dispose(); + gpuContext.Dispose(); + return result; } private static ulong RunHeuristicPipelineCpu( @@ -137,7 +165,7 @@ internal static partial class PrimeOrderCalculator in MontgomeryDivisorData divisorData, ulong phi, PartialFactorResult phiFactors, - PartialFactorGpuContext? gpuContext) + PartialFactorGpuContext gpuContext) { if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) { @@ -388,7 +416,7 @@ internal static partial class PrimeOrderCalculator ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, - PartialFactorGpuContext? gpuContext, + PartialFactorGpuContext gpuContext, out PartialFactorResult? reusableFactorization) { reusableFactorization = null; @@ -551,7 +579,7 @@ internal static partial class PrimeOrderCalculator in PrimeOrderSearchConfig config, PartialFactorResult phiFactors, PartialFactorResult? cachedOrderFactors, - PartialFactorGpuContext? gpuContext, + PartialFactorGpuContext gpuContext, out ulong result) { result = 0UL; @@ -918,7 +946,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, PartialFactorGpuContext? gpuContext, ref int powUsed, int powBudget) + private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, PartialFactorGpuContext gpuContext, ref int powUsed, int powBudget) { PartialFactorResult factorization = PartialFactor(candidate, config, gpuContext); try @@ -1084,7 +1112,7 @@ internal static partial class PrimeOrderCalculator // private static ulong _partialFactorPendingHits; // private static ulong _partialFactorCofactorHits; - private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config, PartialFactorGpuContext? gpuContext = null) + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config, PartialFactorGpuContext gpuContext) { if (value <= 1UL) { @@ -1093,17 +1121,6 @@ internal static partial class PrimeOrderCalculator const int FactorSlotCount = GpuSmallPrimeFactorSlots; - PartialFactorGpuContext? localContext = gpuContext; - bool disposeContext = false; - - if (IsGpuHeuristicDevice && localContext is null) - { - localContext = PartialFactorGpuContext.Rent(); - disposeContext = true; - } - - try - { // stackalloc is faster than pooling // ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); // int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); @@ -1122,30 +1139,21 @@ internal static partial class PrimeOrderCalculator uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; ulong remaining = value; - bool gpuFactored = false; - if (IsGpuHeuristicDevice && localContext is not null) - { - gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor( - localContext.Accelerator, - localContext.Stream, - localContext.Scratch, - value, - limit, - primeSlots, - exponentSlots, - out factorCount, - out remaining, - out _); - } - else if (IsGpuHeuristicDevice) - { - gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); - } + bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor( + gpuContext.Accelerator, + gpuContext.Stream, + gpuContext.Scratch, + value, + limit, + primeSlots, + exponentSlots, + out factorCount, + out remaining, + out _); List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); Stack? compositeStack = null; PartialFactorResult result; - PrimeTester primeTester = PrimeTester.Exclusive; if (!gpuFactored) { @@ -1218,7 +1226,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorHits, 1UL); // Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}"); - bool isPrime = PrimeTester.IsPrimeGpu(composite); + bool isPrime = gpuContext.CheckPrime(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1285,7 +1293,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorPendingHits, 1UL); // Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}"); - bool isPrime = PrimeTester.IsPrimeGpu(composite); + bool isPrime = gpuContext.CheckPrime(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1322,7 +1330,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorCofactorHits, 1UL); // Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}"); - cofactorIsPrime = PrimeTester.IsPrimeCpu(cofactor, CancellationToken.None); + cofactorIsPrime = gpuContext.CheckPrime(cofactor); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1414,14 +1422,6 @@ internal static partial class PrimeOrderCalculator ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); return result; } - finally - { - if (disposeContext && localContext is not null) - { - localContext.Dispose(); - } - } - } private static bool TryPopulateSmallPrimeFactorsCpu( ulong value, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 8e6caf8..89fd553 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -26,50 +26,45 @@ internal static partial class PrimeOrderCalculator { var lease = GpuKernelPool.GetKernel(useGpuOrder: true); var execution = lease.EnterExecutionScope(); - try + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); + + bool success = PrimeOrderGpuHeuristics.TryPartialFactor( + accelerator, + stream, + scratch, + value, + limit, + primeBuffer, + exponentBuffer, + out factorCount, + out remaining, + out _); + + execution.Dispose(); + lease.Dispose(); + + if (!success) { - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); - - bool success = PrimeOrderGpuHeuristics.TryPartialFactor( - accelerator, - stream, - scratch, - value, - limit, - primeBuffer, - exponentBuffer, - out factorCount, - out remaining, - out _); - - if (!success) - { - return false; - } - - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - // This will never happen in production code - // if (primeValue == 0UL || exponent == 0) - // { - // continue; - // } - - counts.Add(primeValue, exponent); - } - - return true; + return false; } - finally + + for (int i = 0; i < factorCount; i++) { - execution.Dispose(); - lease.Dispose(); + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + // This will never happen in production code + // if (primeValue == 0UL || exponent == 0) + // { + // continue; + // } + + counts.Add(primeValue, exponent); } + + return true; } finally { diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index ba46038..2d95f65 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -147,56 +147,12 @@ public sealed class PrimeTester MemoryBuffer1D deviceOutput, Span stagingBuffer) { - if (accelerator is null) - { - throw new ArgumentNullException(nameof(accelerator)); - } - - if (deviceInput is null) - { - throw new ArgumentNullException(nameof(deviceInput)); - } - - if (deviceOutput is null) - { - throw new ArgumentNullException(nameof(deviceOutput)); - } - - if (results.Length < values.Length) - { - throw new ArgumentException("The results span must be at least as long as the input values.", nameof(results)); - } - if (values.IsEmpty) { return; } - if (stagingBuffer.IsEmpty) - { - throw new ArgumentException("The staging buffer must contain at least one element.", nameof(stagingBuffer)); - } - - long inputLength = deviceInput.Length; - long outputLength = deviceOutput.Length; - if (inputLength <= 0) - { - throw new ArgumentException("The device input buffer must contain at least one element.", nameof(deviceInput)); - } - - if (outputLength <= 0) - { - throw new ArgumentException("The device output buffer must contain at least one element.", nameof(deviceOutput)); - } - - long stagingLength = stagingBuffer.Length; - long capacityLong = Math.Min(Math.Min(inputLength, outputLength), stagingLength); - if (capacityLong <= 0) - { - throw new ArgumentException("The provided buffers must be able to hold at least one element."); - } - - int capacity = (int)Math.Min(capacityLong, int.MaxValue); + int capacity = (int)Math.Min(Math.Min(deviceInput.Length, deviceOutput.Length), stagingBuffer.Length); var state = GpuKernelState.GetOrCreate(accelerator); var inputView = deviceInput.View; From 3d458fc6dae28541c92fb157edd7b0ecd1d4748b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 22:28:04 +0100 Subject: [PATCH 4/7] Use precomputed small prime squares in GPU sieve --- .../Gpu/GpuKernelCompilationTests.cs | 4 ++-- .../Gpu/Kernels/PrimeTesterKernels.cs | 9 +++++++-- PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs | 6 ++++-- PerfectNumbers.Core/PrimeTester.cs | 12 +++++++++--- 4 files changed, 22 insertions(+), 9 deletions(-) diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index 1a8b1ae..9f83071 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -121,7 +121,7 @@ public class GpuKernelCompilationTests private static readonly Action, ArrayView1D> MersenneDivisorCyclesKernel = DivisorCycleKernels.GpuDivisorCycleKernel; - private static readonly Action, ArrayView, ArrayView> PrimeTesterSmallPrimeKernel = + private static readonly Action, ArrayView, ArrayView, ArrayView> PrimeTesterSmallPrimeKernel = PrimeTesterKernels.SmallPrimeSieveKernel; private static readonly Action, ArrayView> PrimeTesterSharesFactorKernel = @@ -336,7 +336,7 @@ public class GpuKernelCompilationTests private static void CompilePrimeTesterSmallPrimeKernel(Accelerator accelerator) { - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(PrimeTesterSmallPrimeKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(PrimeTesterSmallPrimeKernel); _ = kernel; } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 1d284e6..9ad7ee7 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -7,7 +7,12 @@ namespace PerfectNumbers.Core.Gpu; internal static class PrimeTesterKernels { // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise. - public static void SmallPrimeSieveKernel(Index1D index, ArrayView numbers, ArrayView smallPrimes, ArrayView results) + public static void SmallPrimeSieveKernel( + Index1D index, + ArrayView numbers, + ArrayView smallPrimes, + ArrayView smallPrimeSquares, + ArrayView results) { ulong n = numbers[index]; @@ -28,7 +33,7 @@ internal static class PrimeTesterKernels for (int i = 0; i < length; i++) { ulong prime = smallPrimes[i]; - ulong primeSquare = (ulong)prime * prime; + ulong primeSquare = smallPrimeSquares[i]; bool withinRange = primeSquare <= n; ulong remainder = n % prime; bool divides = remainder == 0UL; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 4e21cf3..2bbe8b9 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -52,8 +52,9 @@ internal static partial class PrimeOrderCalculator private readonly PrimeTester.KernelState.ScratchBuffers _primeScratch; private readonly ulong[] _primeHostBuffer; private readonly byte[] _primeFlagBuffer; - private readonly Action, ArrayView, ArrayView> _primeKernel; + private readonly Action, ArrayView, ArrayView, ArrayView> _primeKernel; private readonly ArrayView1D _primeDevicePrimes; + private readonly ArrayView1D _primeDevicePrimeSquares; private PartialFactorGpuContext( GpuKernelLease lease, @@ -73,6 +74,7 @@ internal static partial class PrimeOrderCalculator _primeFlagBuffer = primeFlagBuffer; _primeKernel = primeState.Kernel; _primeDevicePrimes = primeState.DevicePrimes.View; + _primeDevicePrimeSquares = primeState.DevicePrimeSquares.View; } public Accelerator Accelerator => _lease.Accelerator; @@ -109,7 +111,7 @@ internal static partial class PrimeOrderCalculator ArrayView1D inputSlice = PrimeDeviceInput.View.SubView(0, 1); inputSlice.CopyFromCPU(ref host[0], 1); ArrayView1D outputSlice = PrimeDeviceOutput.View.SubView(0, 1); - _primeKernel(1, inputSlice.AsContiguous(), _primeDevicePrimes.AsContiguous(), outputSlice.AsContiguous()); + _primeKernel(1, inputSlice.AsContiguous(), _primeDevicePrimes.AsContiguous(), _primeDevicePrimeSquares.AsContiguous(), outputSlice.AsContiguous()); Accelerator.Synchronize(); Span flags = PrimeFlagSpan; outputSlice.CopyToCPU(ref flags[0], 1); diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 2d95f65..09b12c7 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -158,6 +158,7 @@ public sealed class PrimeTester var inputView = deviceInput.View; var outputView = deviceOutput.View; var primesView = state.DevicePrimes.View; + var primeSquaresView = state.DevicePrimeSquares.View; var stagingSpan = stagingBuffer; ref ulong stagingRef = ref MemoryMarshal.GetReference(stagingSpan); @@ -175,7 +176,7 @@ public sealed class PrimeTester inputSlice.CopyFromCPU(ref stagingRef, count); var outputSlice = outputView.SubView(0, count); - state.Kernel(count, inputSlice, primesView, outputSlice); + state.Kernel(count, inputSlice, primesView, primeSquaresView, outputSlice); accelerator.Synchronize(); outputSlice.CopyToCPU(ref results[pos], count); @@ -272,9 +273,10 @@ public sealed class PrimeTester // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). internal sealed class KernelState { - public Action, ArrayView, ArrayView> Kernel { get; } + public Action, ArrayView, ArrayView, ArrayView> Kernel { get; } public Action, ulong, ArrayView> HeuristicTrialDivisionKernel { get; } public MemoryBuffer1D DevicePrimes { get; } + public MemoryBuffer1D DevicePrimeSquares { get; } private readonly Accelerator _accel; private readonly System.Collections.Concurrent.ConcurrentBag _scratchPool = []; // TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in @@ -285,12 +287,15 @@ public sealed class PrimeTester { _accel = accelerator; // Compile once per accelerator and upload primes once. - Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); + Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel, ulong, ArrayView>(PrimeTesterKernels.HeuristicTrialDivisionKernel); var primes = PrimesGenerator.SmallPrimes; + var primeSquares = PrimesGenerator.SmallPrimesPow2; DevicePrimes = accelerator.Allocate1D(primes.Length); + DevicePrimeSquares = accelerator.Allocate1D(primeSquares.Length); DevicePrimes.View.CopyFromCPU(primes); + DevicePrimeSquares.View.CopyFromCPU(primeSquares); } internal sealed class ScratchBuffers : IDisposable @@ -344,6 +349,7 @@ public sealed class PrimeTester } DevicePrimes.Dispose(); + DevicePrimeSquares.Dispose(); } } From 0e1212aa1c1c8c95f0acefe6d6523f03a54ab1e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 22:42:21 +0100 Subject: [PATCH 5/7] Split MulHigh extensions into CPU and GPU variants --- .../Mul64Benchmarks.cs | 2 +- .../MulHighBenchmarks.cs | 10 +++++----- .../MulMod64Benchmarks.cs | 2 +- .../ULongExtensionsTests.cs | 14 ++++++++++++-- PerfectNumbers.Core/UInt128Extensions.cs | 2 +- PerfectNumbers.Core/ULongExtensions.Gpu.cs | 1 + PerfectNumbers.Core/ULongExtensions.cs | 17 +++++++++++------ 7 files changed, 32 insertions(+), 16 deletions(-) diff --git a/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs index 4203709..f5c0d4e 100644 --- a/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs @@ -89,7 +89,7 @@ public class Mul64Benchmarks { ulong leftLow = (ulong)left; ulong rightLow = (ulong)right; - return ((UInt128)(leftLow * (ulong)(right >> 64) + leftLow.MulHigh(rightLow)) << 64) | (UInt128)(leftLow * rightLow); + return ((UInt128)(leftLow * (ulong)(right >> 64) + leftLow.MulHighCpu(rightLow)) << 64) | (UInt128)(leftLow * rightLow); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs index 49e2208..852010c 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs @@ -3,6 +3,7 @@ using System.Runtime.CompilerServices; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; namespace EvenPerfectBitScanner.Benchmarks; @@ -36,10 +37,7 @@ public class MulHighBenchmarks [Benchmark(Baseline = true)] public ulong CurrentImplementation() { - // TODO: Replace ULongExtensions.MulHigh with the UInt128-based helper on CPU paths; it was - // 11–14× faster in the latest benchmarks while the GPU-specific variant remains available - // separately. - return ULongExtensions.MulHigh(Input.X, Input.Y); + return ULongExtensions.MulHighCpu(Input.X, Input.Y); } /// @@ -107,7 +105,9 @@ public class MulHighBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong MulHighWithGpuUInt128(ulong x, ulong y) { - return ULongExtensions.MulHighGpuCompatible(x, y); + GpuUInt128 product = new(x); + GpuUInt128.Mul64(ref product, 0UL, y); + return product.High; } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index d4454a9..949b10a 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -240,7 +240,7 @@ public class MulMod64Benchmarks ulong left = a % modulus; ulong right = b % modulus; ulong low = left * right; - ulong high = left.MulHigh(right); + ulong high = left.MulHighCpu(right); UInt128 product = ((UInt128)high << 64) | low; return (ulong)(product % modulus); } diff --git a/PerfectNumbers.Core.Tests/ULongExtensionsTests.cs b/PerfectNumbers.Core.Tests/ULongExtensionsTests.cs index b192dab..1cbee47 100644 --- a/PerfectNumbers.Core.Tests/ULongExtensionsTests.cs +++ b/PerfectNumbers.Core.Tests/ULongExtensionsTests.cs @@ -131,11 +131,21 @@ public class ULongExtensionsTests [InlineData(0UL, 0UL)] [InlineData(ulong.MaxValue, 1UL)] [InlineData(123456789UL, 987654321UL)] - public void MulHigh_matches_high_bits(ulong x, ulong y) + public void MulHighCpu_matches_high_bits(ulong x, ulong y) { ulong expected = (ulong)(((UInt128)x * y) >> 64); - x.MulHigh(y).Should().Be(expected); + x.MulHighCpu(y).Should().Be(expected); + } + + [Theory] + [Trait("Category", "Fast")] + [InlineData(0UL, 0UL)] + [InlineData(ulong.MaxValue, 1UL)] + [InlineData(123456789UL, 987654321UL)] + public void MulHighGpu_matches_cpu_baseline(ulong x, ulong y) + { + x.MulHighGpu(y).Should().Be(x.MulHighCpu(y)); } [Theory] diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 2e72d58..973bec6 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -760,7 +760,7 @@ public static class UInt128Extensions // expression tree around the shift. Mirroring the MulHigh layout lets RyuJIT // keep the intermediate sum in registers instead of reloading the partial // product from the stack. - ulong high = aLow.MulHigh(bLow); + ulong high = aLow.MulHighCpu(bLow); high += aLow * bHigh; return ((UInt128)high << 64) | low; diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index 6b736ee..2e37a22 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -11,6 +11,7 @@ namespace PerfectNumbers.Core; public static partial class ULongExtensions { [MethodImpl(MethodImplOptions.AggressiveInlining)] + [Obsolete("Use ULongExtensions.MulHighCpu or ULongExtensions.MulHighGpu.")] public static ulong MulHighGpuCompatible(this ulong x, ulong y) { GpuUInt128 product = new(x); diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 08af224..7c3667c 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -430,14 +430,19 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b); + public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHighCpu(b) << 64) | (UInt128)(a * b); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulHigh(this ulong x, ulong y) + public static ulong MulHighCpu(this ulong x, ulong y) { - // TODO: Investigate replacing this manual decomposition with the UInt128-based implementation - // for CPU callers; the latest benchmarks show the intrinsic path is an order of magnitude - // faster, while GPU code can keep using GpuUInt128.MulHigh. + return (ulong)(((UInt128)x * y) >> 64); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulHighGpu(this ulong x, ulong y) + { + // The CPU helpers rely on UInt128 via MulHighCpu; retain this manual decomposition for GPU-style + // arithmetic so callers can avoid the intrinsic when staging work for accelerators. ulong xLow = (uint)x; ulong xHigh = x >> 32; ulong yLow = (uint)y; @@ -554,7 +559,7 @@ public static partial class ULongExtensions ulong m = unchecked(tLow * nPrime); ulong mTimesModulusLow = unchecked(m * modulus); - ulong result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL)); + ulong result = unchecked(a.MulHighCpu(b) + m.MulHighCpu(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL)); if (result >= modulus) { result -= modulus; From e2f68cedf2c2defe22458f184ffb7f8e6d04140b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 22:42:27 +0100 Subject: [PATCH 6/7] Short-circuit small prime sieve kernel loops --- .../Gpu/Kernels/PrimeTesterKernels.cs | 25 ++++++++----------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 9ad7ee7..52a1770 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -15,35 +15,32 @@ internal static class PrimeTesterKernels ArrayView results) { ulong n = numbers[index]; + byte result = 1; // EvenPerfectBitScanner filters candidates so production GPU launches only see: // - n >= 31 (the scanner starts at 136,279,841 and test mode seeds 31) // - n odd (candidate generators and ModResidueTracker eliminate even exponents) // - n not divisible by 5 (residue prefiltering and AddPrimes transform cover these composites) // Tests should honor the same constraints so the host wrapper can remain branchless. - ulong log2Estimate = 63UL - (ulong)XMath.LeadingZeroCount(n | 1UL); - ulong gcd = BinaryGcdGpu(n, log2Estimate); - bool heuristicActive = n.Mod10() == 1UL; - bool gcdReject = heuristicActive & (gcd != 1UL); - - bool composite = gcdReject; - bool continueLoop = !gcdReject; long length = smallPrimes.Length; for (int i = 0; i < length; i++) { ulong prime = smallPrimes[i]; ulong primeSquare = smallPrimeSquares[i]; - bool withinRange = primeSquare <= n; - ulong remainder = n % prime; - bool divides = remainder == 0UL; + if (primeSquare > n) + { + break; + } - bool markComposite = continueLoop & withinRange & divides; - composite |= markComposite; - continueLoop &= withinRange & !divides; + if (n % prime == 0UL) + { + result = 0; + break; + } } - results[index] = composite ? (byte)0 : (byte)1; + results[index] = result; } public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView divisors, ulong n, ArrayView results) From 9a0e32899534636e80a92a18db1b6546b8f832a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 23:41:14 +0100 Subject: [PATCH 7/7] Reintroduce MulHigh CPU helper and tighten GPU sieve loop --- .../MulHighBenchmarks.cs | 6 +- .../Gpu/GpuKernelCompilationTests.cs | 4 +- .../Gpu/PrimeTesterGpuTests.cs | 30 --- PerfectNumbers.Core/Gpu/GpuKernelPool.cs | 16 +- .../Gpu/Kernels/PrimeTesterKernels.cs | 9 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 180 +++++++++++------- .../PrimeOrderCalculator.Cpu.cs | 148 +++----------- .../PrimeOrderCalculator.Gpu.cs | 89 +++++---- PerfectNumbers.Core/PrimeTester.cs | 80 ++------ PerfectNumbers.Core/ULongExtensions.cs | 4 +- 10 files changed, 222 insertions(+), 344 deletions(-) diff --git a/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs index 852010c..a332f23 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs @@ -3,7 +3,6 @@ using System.Runtime.CompilerServices; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using PerfectNumbers.Core; -using PerfectNumbers.Core.Gpu; namespace EvenPerfectBitScanner.Benchmarks; @@ -37,6 +36,7 @@ public class MulHighBenchmarks [Benchmark(Baseline = true)] public ulong CurrentImplementation() { + // Routes through the UInt128-based helper adopted by CPU hot paths after benchmarking. return ULongExtensions.MulHighCpu(Input.X, Input.Y); } @@ -105,9 +105,7 @@ public class MulHighBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong MulHighWithGpuUInt128(ulong x, ulong y) { - GpuUInt128 product = new(x); - GpuUInt128.Mul64(ref product, 0UL, y); - return product.High; + return ULongExtensions.MulHighGpu(x, y); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index 9f83071..1a8b1ae 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -121,7 +121,7 @@ public class GpuKernelCompilationTests private static readonly Action, ArrayView1D> MersenneDivisorCyclesKernel = DivisorCycleKernels.GpuDivisorCycleKernel; - private static readonly Action, ArrayView, ArrayView, ArrayView> PrimeTesterSmallPrimeKernel = + private static readonly Action, ArrayView, ArrayView> PrimeTesterSmallPrimeKernel = PrimeTesterKernels.SmallPrimeSieveKernel; private static readonly Action, ArrayView> PrimeTesterSharesFactorKernel = @@ -336,7 +336,7 @@ public class GpuKernelCompilationTests private static void CompilePrimeTesterSmallPrimeKernel(Accelerator accelerator) { - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(PrimeTesterSmallPrimeKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(PrimeTesterSmallPrimeKernel); _ = kernel; } diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs index 0cce423..5e8faf9 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs @@ -1,6 +1,4 @@ using FluentAssertions; -using ILGPU; -using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; using Xunit; @@ -29,34 +27,6 @@ public class PrimeTesterGpuTests } } - [Fact] - public void IsPrimeBatchGpu_with_preallocated_resources_marks_expected_values() - { - GpuContextPool.ForceCpu = false; - using var context = Context.CreateDefault(); - using var accelerator = context.GetPreferredDevice(false).CreateAccelerator(context); - - try - { - ulong[] values = [31UL, 33UL, 37UL, 39UL]; - byte[] results = new byte[values.Length]; - - using var inputBuffer = accelerator.Allocate1D(values.Length); - using var outputBuffer = accelerator.Allocate1D(values.Length); - ulong[] staging = new ulong[values.Length]; - - PrimeTester.IsPrimeBatchGpu(values, results, accelerator, inputBuffer, outputBuffer, staging); - - results.Should().Equal(new byte[] { 1, 0, 1, 0 }); - } - finally - { - PrimeTester.ClearGpuCaches(accelerator); - GpuContextPool.ForceCpu = false; - GpuContextPool.DisposeAll(); - } - } - [Fact] public void IsPrimeGpu_rejects_composites() { diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 2060b4f..c16d9ca 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -46,7 +46,6 @@ public sealed class KernelContainer public MemoryBuffer1D? SmallPrimeFactorExponentSlots; public MemoryBuffer1D? SmallPrimeFactorCountSlot; public MemoryBuffer1D? SmallPrimeFactorRemainingSlot; - public MemoryBuffer1D? SmallPrimeFactorFullyFactoredSlot; public MemoryBuffer1D? SmallPrimesLastOne; public MemoryBuffer1D? SmallPrimesPow2LastOne; public MemoryBuffer1D? SmallPrimesLastSeven; @@ -108,20 +107,17 @@ public readonly struct SmallPrimeFactorScratch public readonly MemoryBuffer1D ExponentSlots; public readonly MemoryBuffer1D CountSlot; public readonly MemoryBuffer1D RemainingSlot; - public readonly MemoryBuffer1D FullyFactoredSlot; public SmallPrimeFactorScratch( MemoryBuffer1D primeSlots, MemoryBuffer1D exponentSlots, MemoryBuffer1D countSlot, - MemoryBuffer1D remainingSlot, - MemoryBuffer1D fullyFactoredSlot) + MemoryBuffer1D remainingSlot) { PrimeSlots = primeSlots; ExponentSlots = exponentSlots; CountSlot = countSlot; RemainingSlot = remainingSlot; - FullyFactoredSlot = fullyFactoredSlot; } public void Clear() @@ -130,7 +126,6 @@ public readonly struct SmallPrimeFactorScratch ExponentSlots.MemSetToZero(); CountSlot.MemSetToZero(); RemainingSlot.MemSetToZero(); - FullyFactoredSlot.MemSetToZero(); } } @@ -278,14 +273,7 @@ public class GpuKernelPool kernels.SmallPrimeFactorRemainingSlot = remainingSlot; } - MemoryBuffer1D? fullyFactoredSlot = kernels.SmallPrimeFactorFullyFactoredSlot; - if (fullyFactoredSlot is null) - { - fullyFactoredSlot = accelerator.Allocate1D(1); - kernels.SmallPrimeFactorFullyFactoredSlot = fullyFactoredSlot; - } - - return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!, fullyFactoredSlot!); + return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!); } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 52a1770..bf1e357 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -7,12 +7,7 @@ namespace PerfectNumbers.Core.Gpu; internal static class PrimeTesterKernels { // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise. - public static void SmallPrimeSieveKernel( - Index1D index, - ArrayView numbers, - ArrayView smallPrimes, - ArrayView smallPrimeSquares, - ArrayView results) + public static void SmallPrimeSieveKernel(Index1D index, ArrayView numbers, ArrayView smallPrimes, ArrayView results) { ulong n = numbers[index]; byte result = 1; @@ -27,7 +22,7 @@ internal static class PrimeTesterKernels for (int i = 0; i < length; i++) { ulong prime = smallPrimes[i]; - ulong primeSquare = smallPrimeSquares[i]; + ulong primeSquare = (ulong)prime * prime; if (primeSquare > n) { break; diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 45ef028..7b17b48 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -130,11 +130,7 @@ internal static partial class PrimeOrderGpuHeuristics }); } - public static bool TryPartialFactor( - Accelerator accelerator, - AcceleratorStream stream, - SmallPrimeFactorScratch scratch, ulong value, uint limit, Span primeTargets, @@ -147,63 +143,28 @@ internal static partial class PrimeOrderGpuHeuristics remaining = value; fullyFactored = false; + if (primeTargets.Length == 0 || exponentTargets.Length == 0) + { + return false; + } + primeTargets.Clear(); exponentTargets.Clear(); - return TryPartialFactorCore( - accelerator, - stream, - scratch, - value, - limit, - primeTargets, - exponentTargets, - out factorCount, - out remaining, - out fullyFactored); - } + if (!TryLaunchPartialFactorKernel( + value, + limit, + primeTargets, + exponentTargets, + out int extracted, + out ulong leftover, + out bool kernelFullyFactored)) + { + return false; + } - private static bool TryPartialFactorCore( - Accelerator accelerator, - AcceleratorStream stream, - SmallPrimeFactorScratch scratch, - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining, - out bool fullyFactored) - { - factorCount = 0; - remaining = value; - fullyFactored = false; - - scratch.Clear(); - - var kernel = GetPartialFactorKernel(accelerator); - SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); - - kernel( - stream, - 1, - cache.Primes!.View, - cache.Squares!.View, - cache.Count, - primeTargets.Length, - value, - limit, - scratch.PrimeSlots.View, - scratch.ExponentSlots.View, - scratch.CountSlot.View, - scratch.RemainingSlot.View, - scratch.FullyFactoredSlot.View); - - stream.Synchronize(); - - scratch.CountSlot.View.CopyToCPU(ref factorCount, 1); int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); - if (factorCount > capacity) + if (extracted > capacity) { factorCount = 0; remaining = value; @@ -211,21 +172,104 @@ internal static partial class PrimeOrderGpuHeuristics return false; } - if (factorCount > 0) - { - scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeTargets), factorCount); - scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentTargets), factorCount); - } - - scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1); - - byte fullyFactoredFlag = 0; - scratch.FullyFactoredSlot.View.CopyToCPU(ref fullyFactoredFlag, 1); - bool kernelFullyFactored = fullyFactoredFlag != 0; - fullyFactored = kernelFullyFactored && remaining == 1UL; + 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); + var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + var kernel = GetPartialFactorKernel(accelerator); + SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); + + // TODO: We should create / reallocate these buffer only once or if the new required length is bigger than capacity. + var factorBuffer = accelerator.Allocate1D(primeTargets.Length); + var exponentBuffer = accelerator.Allocate1D(exponentTargets.Length); + var countBuffer = accelerator.Allocate1D(1); + var remainingBuffer = accelerator.Allocate1D(1); + var fullyFactoredBuffer = accelerator.Allocate1D(1); + + // TODO: There is no need to clear these buffers because the kernel will always assign values within the required bounds. + // 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; + + // TODO: These buffers shouldn't be disposed but rather reused accross call with the assigned accelerator. + factorBuffer.Dispose(); + exponentBuffer.Dispose(); + countBuffer.Dispose(); + remainingBuffer.Dispose(); + fullyFactoredBuffer.Dispose(); + execution.Dispose(); + lease.Dispose(); + return true; + } + catch (CLException ex) + { + Console.WriteLine($"GPU ERROR ({ex.Error}): {ex.Message}"); + factorCount = 0; + remaining = value; + fullyFactored = false; + return false; + } + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + { + Console.WriteLine($"GPU ERROR: {ex.Message}"); + factorCount = 0; + remaining = value; + fullyFactored = false; + return false; + } + } + public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData divisorData) { Span exponents = stackalloc ulong[1]; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 2bbe8b9..9b83a9c 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -3,10 +3,9 @@ using System.Buffers; using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; -using ILGPU; -using ILGPU.Runtime; -using PerfectNumbers.Core.Gpu; using System.Runtime.InteropServices; +using ILGPU; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -44,90 +43,6 @@ internal static partial class PrimeOrderCalculator } } - private sealed class PartialFactorGpuContext : IDisposable - { - private readonly GpuKernelLease _lease; - private readonly GpuKernelLease.ExecutionScope _executionScope; - private readonly PrimeTester.KernelState _primeState; - private readonly PrimeTester.KernelState.ScratchBuffers _primeScratch; - private readonly ulong[] _primeHostBuffer; - private readonly byte[] _primeFlagBuffer; - private readonly Action, ArrayView, ArrayView, ArrayView> _primeKernel; - private readonly ArrayView1D _primeDevicePrimes; - private readonly ArrayView1D _primeDevicePrimeSquares; - - private PartialFactorGpuContext( - GpuKernelLease lease, - GpuKernelLease.ExecutionScope executionScope, - SmallPrimeFactorScratch scratch, - PrimeTester.KernelState primeState, - PrimeTester.KernelState.ScratchBuffers primeScratch, - ulong[] primeHostBuffer, - byte[] primeFlagBuffer) - { - _lease = lease; - _executionScope = executionScope; - Scratch = scratch; - _primeState = primeState; - _primeScratch = primeScratch; - _primeHostBuffer = primeHostBuffer; - _primeFlagBuffer = primeFlagBuffer; - _primeKernel = primeState.Kernel; - _primeDevicePrimes = primeState.DevicePrimes.View; - _primeDevicePrimeSquares = primeState.DevicePrimeSquares.View; - } - - public Accelerator Accelerator => _lease.Accelerator; - - public AcceleratorStream Stream => _lease.Stream; - - public SmallPrimeFactorScratch Scratch { get; } - - public Span PrimeHostSpan => _primeHostBuffer.AsSpan(0, _primeScratch.Capacity); - - public Span PrimeFlagSpan => _primeFlagBuffer.AsSpan(0, _primeScratch.Capacity); - - public MemoryBuffer1D PrimeDeviceInput => _primeScratch.Input; - - public MemoryBuffer1D PrimeDeviceOutput => _primeScratch.Output; - - public static PartialFactorGpuContext Rent() - { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - var execution = lease.EnterExecutionScope(); - Accelerator accelerator = lease.Accelerator; - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); - PrimeTester.KernelState primeState = PrimeTester.GpuKernelState.GetOrCreate(accelerator); - PrimeTester.KernelState.ScratchBuffers primeScratch = primeState.RentScratch(1, accelerator); - ulong[] primeHostBuffer = ThreadStaticPools.UlongPool.Rent(primeScratch.Capacity); - byte[] primeFlagBuffer = ThreadStaticPools.BytePool.Rent(primeScratch.Capacity); - return new PartialFactorGpuContext(lease, execution, scratch, primeState, primeScratch, primeHostBuffer, primeFlagBuffer); - } - - public bool CheckPrime(ulong value) - { - Span host = PrimeHostSpan; - host[0] = value; - ArrayView1D inputSlice = PrimeDeviceInput.View.SubView(0, 1); - inputSlice.CopyFromCPU(ref host[0], 1); - ArrayView1D outputSlice = PrimeDeviceOutput.View.SubView(0, 1); - _primeKernel(1, inputSlice.AsContiguous(), _primeDevicePrimes.AsContiguous(), _primeDevicePrimeSquares.AsContiguous(), outputSlice.AsContiguous()); - Accelerator.Synchronize(); - Span flags = PrimeFlagSpan; - outputSlice.CopyToCPU(ref flags[0], 1); - return flags[0] != 0; - } - - public void Dispose() - { - _primeState.ReturnScratch(_primeScratch); - ThreadStaticPools.UlongPool.Return(_primeHostBuffer, clearArray: false); - ThreadStaticPools.BytePool.Return(_primeFlagBuffer, clearArray: false); - _executionScope.Dispose(); - _lease.Dispose(); - } - } - private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. @@ -143,20 +58,18 @@ internal static partial class PrimeOrderCalculator return gpuOrder; } - PartialFactorGpuContext gpuContext = PartialFactorGpuContext.Rent(); - PartialFactorResult phiFactors = PartialFactor(phi, config, gpuContext); + PartialFactorResult phiFactors = PartialFactor(phi, config); + ulong result; if (phiFactors.Factors is null) { - ulong fallback = CalculateByFactorizationCpu(prime, divisorData); + result = CalculateByFactorizationCpu(prime, divisorData); phiFactors.Dispose(); - gpuContext.Dispose(); - return fallback; + return result; } - ulong result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors, gpuContext); + result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); phiFactors.Dispose(); - gpuContext.Dispose(); return result; } @@ -166,8 +79,7 @@ internal static partial class PrimeOrderCalculator in PrimeOrderSearchConfig config, in MontgomeryDivisorData divisorData, ulong phi, - PartialFactorResult phiFactors, - PartialFactorGpuContext gpuContext) + PartialFactorResult phiFactors) { if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) { @@ -177,7 +89,7 @@ internal static partial class PrimeOrderCalculator ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, gpuContext, out PartialFactorResult? orderFactors)) + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) { return candidateOrder; } @@ -188,7 +100,7 @@ internal static partial class PrimeOrderCalculator return CalculateByFactorizationCpu(prime, divisorData); } - if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, gpuContext, out ulong order)) + if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) { return order; } @@ -418,7 +330,6 @@ internal static partial class PrimeOrderCalculator ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, - PartialFactorGpuContext gpuContext, out PartialFactorResult? reusableFactorization) { reusableFactorization = null; @@ -433,7 +344,7 @@ internal static partial class PrimeOrderCalculator return false; } - PartialFactorResult factorization = PartialFactor(order, config, gpuContext); + PartialFactorResult factorization = PartialFactor(order, config); if (factorization.Factors is null) { factorization.Dispose(); @@ -581,7 +492,6 @@ internal static partial class PrimeOrderCalculator in PrimeOrderSearchConfig config, PartialFactorResult phiFactors, PartialFactorResult? cachedOrderFactors, - PartialFactorGpuContext gpuContext, out ulong result) { result = 0UL; @@ -592,7 +502,7 @@ internal static partial class PrimeOrderCalculator } // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactor(order, config, gpuContext); + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactor(order, config); try { if (orderFactors.Factors is null) @@ -708,7 +618,7 @@ internal static partial class PrimeOrderCalculator bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); powStepperInitialized = true; - if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, gpuContext, ref powUsed, powBudget)) + if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) { if (gpuPool is not null) { @@ -745,7 +655,7 @@ internal static partial class PrimeOrderCalculator continue; } - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, gpuContext, ref powUsed, powBudget)) + if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) { continue; } @@ -948,9 +858,9 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, PartialFactorGpuContext gpuContext, ref int powUsed, int powBudget) + private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) { - PartialFactorResult factorization = PartialFactor(candidate, config, gpuContext); + PartialFactorResult factorization = PartialFactor(candidate, config); try { if (factorization.Factors is null) @@ -1114,7 +1024,7 @@ internal static partial class PrimeOrderCalculator // private static ulong _partialFactorPendingHits; // private static ulong _partialFactorCofactorHits; - private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config, PartialFactorGpuContext gpuContext) + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) { if (value <= 1UL) { @@ -1141,21 +1051,16 @@ internal static partial class PrimeOrderCalculator uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; ulong remaining = value; - bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor( - gpuContext.Accelerator, - gpuContext.Stream, - gpuContext.Scratch, - value, - limit, - primeSlots, - exponentSlots, - out factorCount, - out remaining, - out _); + bool gpuFactored = false; + if (IsGpuHeuristicDevice) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); + } List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); Stack? compositeStack = null; PartialFactorResult result; + PrimeTester primeTester = PrimeTester.Exclusive; if (!gpuFactored) { @@ -1228,7 +1133,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorHits, 1UL); // Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}"); - bool isPrime = gpuContext.CheckPrime(composite); + bool isPrime = PrimeTester.IsPrimeGpu(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1295,7 +1200,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorPendingHits, 1UL); // Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}"); - bool isPrime = gpuContext.CheckPrime(composite); + bool isPrime = PrimeTester.IsPrimeGpu(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1332,7 +1237,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorCofactorHits, 1UL); // Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}"); - cofactorIsPrime = gpuContext.CheckPrime(cofactor); + cofactorIsPrime = PrimeTester.IsPrimeCpu(cofactor, CancellationToken.None); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1425,6 +1330,7 @@ internal static partial class PrimeOrderCalculator return result; } + private static bool TryPopulateSmallPrimeFactorsCpu( ulong value, uint limit, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 89fd553..4cb801c 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,4 +1,5 @@ using System; +using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; using ILGPU.Runtime.OpenCL; @@ -26,45 +27,63 @@ internal static partial class PrimeOrderCalculator { var lease = GpuKernelPool.GetKernel(useGpuOrder: true); var execution = lease.EnterExecutionScope(); - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - - SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots); - - bool success = PrimeOrderGpuHeuristics.TryPartialFactor( - accelerator, - stream, - scratch, - value, - limit, - primeBuffer, - exponentBuffer, - out factorCount, - out remaining, - out _); - - execution.Dispose(); - lease.Dispose(); - - if (!success) + try { - return false; - } + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; - for (int i = 0; i < factorCount; i++) + 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(); + + 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]; + // This will never happen in production code + // if (primeValue == 0UL || exponent == 0) + // { + // continue; + // } + + counts.Add(primeValue, exponent); + } + + return true; + } + finally { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - // This will never happen in production code - // if (primeValue == 0UL || exponent == 0) - // { - // continue; - // } - - counts.Add(primeValue, exponent); + execution.Dispose(); + lease.Dispose(); } - - return true; } finally { diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 09b12c7..dd15c4e 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -3,7 +3,6 @@ using System.Buffers; using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; -using System.Runtime.InteropServices; using System.Threading; using ILGPU; using ILGPU.Algorithms; @@ -113,6 +112,7 @@ public sealed class PrimeTester var gpu = PrimeTesterGpuContextPool.Rent(); var accelerator = gpu.Accelerator; var state = GpuKernelState.GetOrCreate(accelerator); + int totalLength = values.Length; int batchSize = Math.Max(1, GpuBatchSize); lock (gpu.ExecutionLock) @@ -123,67 +123,30 @@ public sealed class PrimeTester ArrayPool pool = ThreadStaticPools.UlongPool; ulong[] temp = pool.Rent(batchSize); - try + int pos = 0; + while (pos < totalLength) { - IsPrimeBatchGpu(values, results, accelerator, input, output, temp); - } - finally - { - pool.Return(temp, clearArray: false); - state.ReturnScratch(scratch); + int remaining = totalLength - pos; + int count = remaining > batchSize ? batchSize : remaining; + + values.Slice(pos, count).CopyTo(temp); + input.View.CopyFromCPU(ref temp[0], count); + + state.Kernel(count, input.View, state.DevicePrimes.View, output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref results[pos], count); + + pos += count; } + + pool.Return(temp, clearArray: false); + state.ReturnScratch(scratch); } gpu.Dispose(); limiter.Dispose(); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static void IsPrimeBatchGpu( - ReadOnlySpan values, - Span results, - Accelerator accelerator, - MemoryBuffer1D deviceInput, - MemoryBuffer1D deviceOutput, - Span stagingBuffer) - { - if (values.IsEmpty) - { - return; - } - - int capacity = (int)Math.Min(Math.Min(deviceInput.Length, deviceOutput.Length), stagingBuffer.Length); - - var state = GpuKernelState.GetOrCreate(accelerator); - var inputView = deviceInput.View; - var outputView = deviceOutput.View; - var primesView = state.DevicePrimes.View; - var primeSquaresView = state.DevicePrimeSquares.View; - var stagingSpan = stagingBuffer; - ref ulong stagingRef = ref MemoryMarshal.GetReference(stagingSpan); - - int totalLength = values.Length; - int pos = 0; - while (pos < totalLength) - { - int remaining = totalLength - pos; - int count = remaining > capacity ? capacity : remaining; - - var stagingSlice = stagingSpan.Slice(0, count); - values.Slice(pos, count).CopyTo(stagingSlice); - - var inputSlice = inputView.SubView(0, count); - inputSlice.CopyFromCPU(ref stagingRef, count); - - var outputSlice = outputView.SubView(0, count); - state.Kernel(count, inputSlice, primesView, primeSquaresView, outputSlice); - accelerator.Synchronize(); - outputSlice.CopyToCPU(ref results[pos], count); - - pos += count; - } - } - internal static class PrimeTesterGpuContextPool { internal sealed class PooledContext : IDisposable @@ -273,10 +236,9 @@ public sealed class PrimeTester // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). internal sealed class KernelState { - public Action, ArrayView, ArrayView, ArrayView> Kernel { get; } + public Action, ArrayView, ArrayView> Kernel { get; } public Action, ulong, ArrayView> HeuristicTrialDivisionKernel { get; } public MemoryBuffer1D DevicePrimes { get; } - public MemoryBuffer1D DevicePrimeSquares { get; } private readonly Accelerator _accel; private readonly System.Collections.Concurrent.ConcurrentBag _scratchPool = []; // TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in @@ -287,15 +249,12 @@ public sealed class PrimeTester { _accel = accelerator; // Compile once per accelerator and upload primes once. - Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); + Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel, ulong, ArrayView>(PrimeTesterKernels.HeuristicTrialDivisionKernel); var primes = PrimesGenerator.SmallPrimes; - var primeSquares = PrimesGenerator.SmallPrimesPow2; DevicePrimes = accelerator.Allocate1D(primes.Length); - DevicePrimeSquares = accelerator.Allocate1D(primeSquares.Length); DevicePrimes.View.CopyFromCPU(primes); - DevicePrimeSquares.View.CopyFromCPU(primeSquares); } internal sealed class ScratchBuffers : IDisposable @@ -349,7 +308,6 @@ public sealed class PrimeTester } DevicePrimes.Dispose(); - DevicePrimeSquares.Dispose(); } } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 7c3667c..4c057ef 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -441,8 +441,8 @@ public static partial class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulHighGpu(this ulong x, ulong y) { - // The CPU helpers rely on UInt128 via MulHighCpu; retain this manual decomposition for GPU-style - // arithmetic so callers can avoid the intrinsic when staging work for accelerators. + // Retain this manual decomposition for GPU-style arithmetic so callers staging accelerator work + // can avoid the UInt128 intrinsic while keeping parity with GpuUInt128.MulHigh. ulong xLow = (uint)x; ulong xHigh = x >> 32; ulong yLow = (uint)y;