diff --git a/Benchmarks.txt b/Benchmarks.txt index f26ba28..7cd458d 100644 --- a/Benchmarks.txt +++ b/Benchmarks.txt @@ -23,7 +23,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 |------------------------------------------- |------------------ |-----------:|----------:|----------:|------:|--------:|----------:|------------:| | ExtensionBaseline | CrossWordBlend | 3.606 ns | 0.0030 ns | 0.0025 ns | 1.00 | 0.00 | - | NA | | GpuCompatibleMulModExtension | CrossWordBlend | 22.870 ns | 0.1463 ns | 0.1297 ns | 6.34 | 0.04 | - | NA | -| GpuCompatibleMulModSimplifiedExtension | CrossWordBlend | 22.780 ns | 0.4463 ns | 0.8600 ns | 6.32 | 0.23 | - | NA | +| MulModGpu | CrossWordBlend | 22.780 ns | 0.4463 ns | 0.8600 ns | 6.32 | 0.23 | - | NA | | InlineUInt128Operands | CrossWordBlend | 22.218 ns | 0.0319 ns | 0.0266 ns | 6.16 | 0.01 | - | NA | | InlineUInt128OperandsWithReductionFirst | CrossWordBlend | 22.626 ns | 0.0553 ns | 0.0517 ns | 6.27 | 0.02 | - | NA | | InlineUInt128OperandsWithOperandReduction | CrossWordBlend | 22.548 ns | 0.0487 ns | 0.0432 ns | 6.25 | 0.01 | - | NA | @@ -35,7 +35,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 | | | | | | | | | | | ExtensionBaseline | MixedBitPattern | 3.564 ns | 0.0059 ns | 0.0052 ns | 1.00 | 0.00 | - | NA | | GpuCompatibleMulModExtension | MixedBitPattern | 64.800 ns | 0.3399 ns | 0.3179 ns | 18.17 | 0.09 | - | NA | -| GpuCompatibleMulModSimplifiedExtension | MixedBitPattern | 59.165 ns | 0.4243 ns | 0.3969 ns | 16.60 | 0.11 | - | NA | +| MulModGpu | MixedBitPattern | 59.165 ns | 0.4243 ns | 0.3969 ns | 16.60 | 0.11 | - | NA | | InlineUInt128Operands | MixedBitPattern | 30.539 ns | 0.1396 ns | 0.1306 ns | 8.57 | 0.03 | - | NA | | InlineUInt128OperandsWithReductionFirst | MixedBitPattern | 31.237 ns | 0.1013 ns | 0.0898 ns | 8.76 | 0.03 | - | NA | | InlineUInt128OperandsWithOperandReduction | MixedBitPattern | 30.210 ns | 0.1323 ns | 0.1237 ns | 8.48 | 0.04 | - | NA | @@ -47,7 +47,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 | | | | | | | | | | | ExtensionBaseline | NearFullRange | 3.583 ns | 0.0105 ns | 0.0098 ns | 1.00 | 0.00 | - | NA | | GpuCompatibleMulModExtension | NearFullRange | 5.708 ns | 0.0520 ns | 0.0486 ns | 1.59 | 0.02 | - | NA | -| GpuCompatibleMulModSimplifiedExtension | NearFullRange | 5.952 ns | 0.1253 ns | 0.1172 ns | 1.66 | 0.04 | - | NA | +| MulModGpu | NearFullRange | 5.952 ns | 0.1253 ns | 0.1172 ns | 1.66 | 0.04 | - | NA | | InlineUInt128Operands | NearFullRange | 27.159 ns | 0.0936 ns | 0.0781 ns | 7.58 | 0.03 | - | NA | | InlineUInt128OperandsWithReductionFirst | NearFullRange | 4.357 ns | 0.0205 ns | 0.0182 ns | 1.22 | 0.00 | - | NA | | InlineUInt128OperandsWithOperandReduction | NearFullRange | 4.319 ns | 0.0158 ns | 0.0140 ns | 1.21 | 0.00 | - | NA | @@ -59,7 +59,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 | | | | | | | | | | | ExtensionBaseline | PrimeSizedModulus | 3.593 ns | 0.0039 ns | 0.0037 ns | 1.00 | 0.00 | - | NA | | GpuCompatibleMulModExtension | PrimeSizedModulus | 52.766 ns | 0.3086 ns | 0.2887 ns | 14.69 | 0.08 | - | NA | -| GpuCompatibleMulModSimplifiedExtension | PrimeSizedModulus | 50.490 ns | 0.1977 ns | 0.1752 ns | 14.05 | 0.05 | - | NA | +| MulModGpu | PrimeSizedModulus | 50.490 ns | 0.1977 ns | 0.1752 ns | 14.05 | 0.05 | - | NA | | InlineUInt128Operands | PrimeSizedModulus | 31.409 ns | 0.0984 ns | 0.0920 ns | 8.74 | 0.03 | - | NA | | InlineUInt128OperandsWithReductionFirst | PrimeSizedModulus | 32.176 ns | 0.3427 ns | 0.3206 ns | 8.96 | 0.09 | - | NA | | InlineUInt128OperandsWithOperandReduction | PrimeSizedModulus | 31.786 ns | 0.1483 ns | 0.1387 ns | 8.85 | 0.04 | - | NA | @@ -71,7 +71,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 | | | | | | | | | | | ExtensionBaseline | SparseOperands | 3.555 ns | 0.0087 ns | 0.0082 ns | 1.00 | 0.00 | - | NA | | GpuCompatibleMulModExtension | SparseOperands | 5.671 ns | 0.0043 ns | 0.0036 ns | 1.60 | 0.00 | - | NA | -| GpuCompatibleMulModSimplifiedExtension | SparseOperands | 5.239 ns | 0.0103 ns | 0.0096 ns | 1.47 | 0.00 | - | NA | +| MulModGpu | SparseOperands | 5.239 ns | 0.0103 ns | 0.0096 ns | 1.47 | 0.00 | - | NA | | InlineUInt128Operands | SparseOperands | 3.563 ns | 0.0217 ns | 0.0203 ns | 1.00 | 0.01 | - | NA | | InlineUInt128OperandsWithReductionFirst | SparseOperands | 4.617 ns | 0.0095 ns | 0.0089 ns | 1.30 | 0.00 | - | NA | | InlineUInt128OperandsWithOperandReduction | SparseOperands | 4.282 ns | 0.0193 ns | 0.0181 ns | 1.20 | 0.01 | - | NA | @@ -83,7 +83,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 | | | | | | | | | | | ExtensionBaseline | ZeroOperands | 3.599 ns | 0.0069 ns | 0.0065 ns | 1.00 | 0.00 | - | NA | | GpuCompatibleMulModExtension | ZeroOperands | 1.976 ns | 0.0053 ns | 0.0049 ns | 0.55 | 0.00 | - | NA | -| GpuCompatibleMulModSimplifiedExtension | ZeroOperands | 1.986 ns | 0.0048 ns | 0.0045 ns | 0.55 | 0.00 | - | NA | +| MulModGpu | ZeroOperands | 1.986 ns | 0.0048 ns | 0.0045 ns | 0.55 | 0.00 | - | NA | | InlineUInt128Operands | ZeroOperands | 3.542 ns | 0.0125 ns | 0.0117 ns | 0.98 | 0.00 | - | NA | | InlineUInt128OperandsWithReductionFirst | ZeroOperands | 4.286 ns | 0.0178 ns | 0.0167 ns | 1.19 | 0.01 | - | NA | | InlineUInt128OperandsWithOperandReduction | ZeroOperands | 4.256 ns | 0.0116 ns | 0.0097 ns | 1.18 | 0.00 | - | NA | @@ -97,7 +97,7 @@ Job=.NET 8.0 Runtime=.NET 8.0 Outliers MulMod64Benchmarks.ExtensionBaseline: .NET 8.0 -> 2 outliers were removed (4.92 ns, 5.01 ns) MulMod64Benchmarks.GpuCompatibleMulModExtension: .NET 8.0 -> 1 outlier was removed, 2 outliers were detected (23.75 ns, 25.37 ns) - MulMod64Benchmarks.GpuCompatibleMulModSimplifiedExtension: .NET 8.0 -> 2 outliers were removed, 3 outliers were detected (22.07 ns, 26.53 ns, 26.65 ns) + MulMod64Benchmarks.MulModGpu: .NET 8.0 -> 2 outliers were removed, 3 outliers were detected (22.07 ns, 26.53 ns, 26.65 ns) MulMod64Benchmarks.InlineUInt128Operands: .NET 8.0 -> 2 outliers were removed (23.62 ns, 23.64 ns) MulMod64Benchmarks.InlineUInt128OperandsWithOperandReduction: .NET 8.0 -> 1 outlier was removed (24.42 ns) MulMod64Benchmarks.MultiplyHighDecomposition: .NET 8.0 -> 2 outliers were removed (25.12 ns, 25.31 ns) @@ -112,7 +112,7 @@ Outliers MulMod64Benchmarks.InlineUInt128Operands: .NET 8.0 -> 2 outliers were removed (29.01 ns, 29.49 ns) MulMod64Benchmarks.InlineUInt128OperandsWithReductionFirst: .NET 8.0 -> 1 outlier was removed (5.72 ns) MulMod64Benchmarks.InlineUInt128OperandsWithOperandReduction: .NET 8.0 -> 1 outlier was removed (5.73 ns) - MulMod64Benchmarks.GpuCompatibleMulModSimplifiedExtension: .NET 8.0 -> 1 outlier was removed, 3 outliers were detected (51.41 ns, 51.47 ns, 52.10 ns) + MulMod64Benchmarks.MulModGpu: .NET 8.0 -> 1 outlier was removed, 3 outliers were detected (51.41 ns, 51.47 ns, 52.10 ns) MulMod64Benchmarks.InlineUInt128Operands: .NET 8.0 -> 1 outlier was detected (32.44 ns) MulMod64Benchmarks.InlineUInt128WithLocals: .NET 8.0 -> 2 outliers were removed (33.47 ns, 34.46 ns) MulMod64Benchmarks.MultiplyHighDecomposition: .NET 8.0 -> 1 outlier was removed (33.93 ns) diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.Gpu.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.Gpu.cs new file mode 100644 index 0000000..a7f5d28 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.Gpu.cs @@ -0,0 +1,461 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Runtime.CompilerServices; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +/// +/// GPU benchmarks for 64-bit modular multiplication helpers. +/// Windows 11 results ranged from 0.62 ms to 1.63 ms, and MulModGpu led the PrimeSizedModulus case at 0.641 ms. +/// Recommended helper for p ≥ 138,000,000 is MulModGpu. +/// +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class MulMod64BenchmarksGpu : IDisposable +{ + private static readonly MulMod64Input[] Inputs = + [ + new(0UL, 0UL, 3UL, "ZeroOperands"), + new(ulong.MaxValue, ulong.MaxValue - 1UL, ulong.MaxValue - 58UL, "NearFullRange"), + new(0xFFFF_FFFF_0000_0000UL, 0x0000_0000_FFFF_FFFFUL, 0xFFFF_FFFF_0000_002BUL, "CrossWordBlend"), + new(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL, 0x1FFF_FFFF_FFFF_FFFBUL, "MixedBitPattern"), + new(0x0000_0001_0000_0000UL, 0x0000_0000_0000_0003UL, 0x0000_0000_FFFF_FFC3UL, "SparseOperands"), + new(0x7FFF_FFFF_FFFF_FFA3UL, 0x6FFF_FFFF_FFFF_FF81UL, 0x7FFF_FFFF_FFFF_FFE7UL, "PrimeSizedModulus"), + ]; + + [ParamsSource(nameof(GetInputs))] + public MulMod64Input Input { get; set; } + + private Context? _context; + private Accelerator? _accelerator; + private MemoryBuffer1D? _inputBuffer; + private MemoryBuffer1D? _resultBuffer; + + private Action, ArrayView>? _extensionBaselineKernel; + private Action, ArrayView>? _gpuCompatibleExtensionKernel; + private Action, ArrayView>? _gpuCompatibleSimplifiedKernel; + private Action, ArrayView>? _inlineKernel; + private Action, ArrayView>? _inlineReductionFirstKernel; + private Action, ArrayView>? _inlineOperandReductionKernel; + private Action, ArrayView>? _inlineLocalsKernel; + private Action, ArrayView>? _inlineLocalsOperandReductionKernel; + private Action, ArrayView>? _mulHighKernel; + private Action, ArrayView>? _legacyKernel; + private Action, ArrayView>? _legacyDeferredKernel; + + public static IEnumerable GetInputs() + { + return Inputs; + } + + [GlobalSetup] + public void GlobalSetup() + { + _context = Context.Create(builder => builder.Default().EnableAlgorithms()); + _accelerator = _context.GetPreferredDevice(preferCPU: false).CreateAccelerator(_context); + + MulMod64GpuKernelInput[] gpuInputs = Inputs + .Select(i => new MulMod64GpuKernelInput(i.Left, i.Right, i.Modulus)) + .ToArray(); + + _inputBuffer = _accelerator.Allocate1D(gpuInputs); + _resultBuffer = _accelerator.Allocate1D(gpuInputs.Length); + + _extensionBaselineKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(ExtensionBaselineKernel); + _gpuCompatibleExtensionKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(GpuCompatibleExtensionKernel); + _gpuCompatibleSimplifiedKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(GpuCompatibleSimplifiedKernel); + _inlineKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(InlineKernel); + _inlineReductionFirstKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(InlineReductionFirstKernel); + _inlineOperandReductionKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(InlineOperandReductionKernel); + _inlineLocalsKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(InlineLocalsKernel); + _inlineLocalsOperandReductionKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(InlineLocalsOperandReductionKernel); + _mulHighKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(MulHighKernel); + _legacyKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(LegacyKernel); + _legacyDeferredKernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(LegacyDeferredKernel); + } + + /// + /// Inline wide multiply/reduce baseline: ZeroOperands 0.618 ms, MixedBitPattern 0.640 ms, CrossWordBlend 0.928 ms, NearFullRange 0.875 ms, PrimeSizedModulus 0.759 ms, SparseOperands 0.804 ms. + /// Use primarily as a reference; optimized GPU-specific helpers beat it across all datasets. + /// + [Benchmark(Baseline = true)] + public void ExtensionBaseline() + { + LaunchKernel(_extensionBaselineKernel); + } + + /// + /// Original GPU-compatible extension path: means from 0.614 ms (MixedBitPattern) to 1.023 ms (CrossWordBlend) with PrimeSizedModulus at 0.697 ms. + /// Retain only when mirroring the CPU extension logic outweighs the speed gap to the simplified helper. + /// + [Benchmark] + public void GpuCompatibleMulModExtension() + { + LaunchKernel(_gpuCompatibleExtensionKernel); + } + + /// + /// Simplified GPU extension helper: 0.708 ms (ZeroOperands), 0.905 ms (CrossWordBlend), 0.798 ms (NearFullRange), and class-leading 0.641 ms (PrimeSizedModulus). + /// Recommended for production GPU scans, especially for large p ≥ 138,000,000. + /// + [Benchmark] + public void MulModGpu() + { + LaunchKernel(_gpuCompatibleSimplifiedKernel); + } + + /// + /// Inline UInt128 multiplication without operand reduction: CrossWordBlend 0.838 ms, MixedBitPattern 0.865 ms, NearFullRange 0.786 ms, PrimeSizedModulus 0.709 ms, SparseOperands 0.751 ms, ZeroOperands 0.827 ms. + /// Suitable when explicit UInt128 math is acceptable, though it trails the simplified GPU helper on every dataset. + /// + [Benchmark] + public void InlineUInt128Operands() + { + LaunchKernel(_inlineKernel); + } + + /// + /// Reduces operands before multiplying: ZeroOperands 0.926 ms, CrossWordBlend 0.847 ms, MixedBitPattern 0.843 ms, NearFullRange 0.689 ms, PrimeSizedModulus 0.809 ms, SparseOperands 0.746 ms. + /// Only competitive on NearFullRange; the extra reductions hurt other scenarios. + /// + [Benchmark] + public void InlineUInt128OperandsWithReductionFirst() + { + LaunchKernel(_inlineReductionFirstKernel); + } + + /// + /// Performs operand reduction and then multiplies inline: CrossWordBlend 1.041 ms, MixedBitPattern 0.913 ms, NearFullRange 0.847 ms, PrimeSizedModulus 0.812 ms, SparseOperands 0.744 ms, ZeroOperands 0.802 ms. + /// Consistently slower than both the reduction-first and simplified GPU helpers. + /// + [Benchmark] + public void InlineUInt128OperandsWithOperandReduction() + { + LaunchKernel(_inlineOperandReductionKernel); + } + + /// + /// Uses local temporaries for the wide product: ZeroOperands 0.691 ms, CrossWordBlend 0.914 ms, MixedBitPattern 0.874 ms, NearFullRange 0.702 ms, PrimeSizedModulus 0.817 ms, SparseOperands 0.728 ms. + /// Balanced choice when sticking with manual wide products, but still slower than the simplified GPU path. + /// + [Benchmark] + public void InlineUInt128WithLocals() + { + LaunchKernel(_inlineLocalsKernel); + } + + /// + /// Combines local temporaries with operand reduction: ZeroOperands 0.698 ms, CrossWordBlend 0.858 ms, MixedBitPattern 0.877 ms, NearFullRange 0.851 ms, PrimeSizedModulus 0.761 ms, SparseOperands 0.779 ms. + /// Adds little benefit over the non-reduced local variant while increasing cost on most datasets. + /// + [Benchmark] + public void InlineUInt128WithLocalsAndOperandReduction() + { + LaunchKernel(_inlineLocalsOperandReductionKernel); + } + + /// + /// Multiply-high decomposition: ZeroOperands 0.651 ms, CrossWordBlend 0.848 ms, MixedBitPattern 0.758 ms, NearFullRange 0.812 ms, PrimeSizedModulus 0.728 ms, SparseOperands 0.699 ms. + /// Useful when avoiding UInt128 temporaries on sparse inputs, but still behind the simplified GPU helper for prime-sized moduli. + /// + [Benchmark] + public void MultiplyHighDecomposition() + { + LaunchKernel(_mulHighKernel); + } + + /// + /// Legacy GPU baseline: ZeroOperands 0.759 ms, CrossWordBlend 0.971 ms, MixedBitPattern 0.827 ms, NearFullRange 0.870 ms, PrimeSizedModulus 0.737 ms, SparseOperands 0.620 ms. + /// Retain only for regression comparisons; newer helpers are faster except on the sparse sample. + /// + [Benchmark] + public void GpuCompatibleBaseline() + { + LaunchKernel(_legacyKernel); + } + + /// + /// Deferred legacy helper: ZeroOperands 1.561 ms, CrossWordBlend 1.631 ms, MixedBitPattern 1.540 ms, NearFullRange 1.438 ms, PrimeSizedModulus 1.598 ms, SparseOperands 1.096 ms. + /// Slowest option in every scenario; keep only for historical context. + /// + [Benchmark] + public void GpuCompatibleDeferred() + { + LaunchKernel(_legacyDeferredKernel); + } + + [GlobalCleanup] + public void GlobalCleanup() + { + _resultBuffer?.Dispose(); + _inputBuffer?.Dispose(); + _accelerator?.Dispose(); + _context?.Dispose(); + } + + public void Dispose() + { + GlobalCleanup(); + GC.SuppressFinalize(this); + } + + private void LaunchKernel(Action, ArrayView>? kernel) + { + if (kernel == null || _accelerator == null || _inputBuffer == null || _resultBuffer == null) + { + throw new InvalidOperationException("Benchmark kernels are not initialized."); + } + + kernel(new Index1D((int)_inputBuffer.Length), _inputBuffer.View, _resultBuffer.View); + _accelerator.Synchronize(); + } + + private static void ExtensionBaselineKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + results[index] = MulMod64Baseline(input.Left, input.Right, input.Modulus); + } + + private static void GpuCompatibleExtensionKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + GpuUInt128 state = new(input.Left); + results[index] = state.MulMod(input.Right, input.Modulus); + } + + private static void GpuCompatibleSimplifiedKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + GpuUInt128 state = new(input.Left); + results[index] = state.MulMod(input.Right, input.Modulus); + } + + private static void InlineKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + results[index] = MulMod64Baseline(input.Left, input.Right, input.Modulus); + } + + private static void InlineReductionFirstKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + ulong left = input.Left % input.Modulus; + ulong right = input.Right % input.Modulus; + results[index] = MulMod64Baseline(left, right, input.Modulus); + } + + private static void InlineOperandReductionKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + ulong left = input.Left % input.Modulus; + ulong right = input.Right % input.Modulus; + results[index] = MulMod64Baseline(left, right, input.Modulus); + } + + private static void InlineLocalsKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + MultiplyWide(input.Left, input.Right, out ulong productHigh, out ulong productLow); + results[index] = ReduceWide(productHigh, productLow, input.Modulus); + } + + private static void InlineLocalsOperandReductionKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + ulong left = input.Left % input.Modulus; + ulong right = input.Right % input.Modulus; + MultiplyWide(left, right, out ulong productHigh, out ulong productLow); + results[index] = ReduceWide(productHigh, productLow, input.Modulus); + } + + private static void MulHighKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + ulong left = input.Left % input.Modulus; + ulong right = input.Right % input.Modulus; + MultiplyWide(left, right, out ulong productHigh, out ulong productLow); + results[index] = ReduceWide(productHigh, productLow, input.Modulus); + } + + private static void LegacyKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + results[index] = MulMod64GpuLegacy(input.Left, input.Right, input.Modulus); + } + + private static void LegacyDeferredKernel(Index1D index, ArrayView inputs, ArrayView results) + { + if (index >= inputs.Length) + { + return; + } + + MulMod64GpuKernelInput input = inputs[index]; + results[index] = MulMod64GpuDeferredLegacy(input.Left, input.Right, input.Modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong MulMod64Baseline(ulong left, ulong right, ulong modulus) + { + MultiplyWide(left, right, out ulong high, out ulong low); + return ReduceWide(high, low, modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void MultiplyWide(ulong left, ulong right, out ulong high, out ulong low) + { + ulong leftLow = (uint)left; + ulong leftHigh = left >> 32; + ulong rightLow = (uint)right; + ulong rightHigh = right >> 32; + + ulong productLL = leftLow * rightLow; + ulong productLH = leftLow * rightHigh; + ulong productHL = leftHigh * rightLow; + ulong productHH = leftHigh * rightHigh; + + ulong carry = (productLL >> 32) + (productLH & 0xFFFF_FFFFUL) + (productHL & 0xFFFF_FFFFUL); + low = (carry << 32) | (uint)productLL; + high = productHH + (productLH >> 32) + (productHL >> 32) + (carry >> 32); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ReduceWide(ulong high, ulong low, ulong modulus) + { + if (modulus <= 1UL) + { + return 0UL; + } + + ulong remainder = 0UL; + if (high != 0UL) + { + for (int bit = 63; bit >= 0; --bit) + { + remainder = ShiftLeftAddBit(remainder, (high >> bit) & 1UL, modulus); + } + } + + for (int bit = 63; bit >= 0; --bit) + { + remainder = ShiftLeftAddBit(remainder, (low >> bit) & 1UL, modulus); + } + + return remainder; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ShiftLeftAddBit(ulong value, ulong bit, ulong modulus) + { + value <<= 1; + if (value >= modulus) + { + value -= modulus; + } + + if (bit != 0UL) + { + value += 1UL; + if (value >= modulus) + { + value -= modulus; + } + } + + return value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong MulMod64GpuLegacy(ulong a, ulong b, ulong modulus) + { + GpuUInt128 state = new(a % modulus); + return state.MulMod(b, modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong MulMod64GpuDeferredLegacy(ulong a, ulong b, ulong modulus) + { + GpuUInt128 state = new(a); + return state.MulModWithNativeModulo(b, modulus); + } + + public readonly record struct MulMod64Input(ulong Left, ulong Right, ulong Modulus, string Name) + { + public override string ToString() + { + return Name; + } + } + + private readonly struct MulMod64GpuKernelInput + { + public readonly ulong Left; + public readonly ulong Right; + public readonly ulong Modulus; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public MulMod64GpuKernelInput(ulong left, ulong right, ulong modulus) + { + Left = left; + Right = right; + Modulus = modulus; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index beb7798..d4454a9 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -58,7 +58,7 @@ public class MulMod64Benchmarks } /// - /// Calls the UInt128-based helper; performance mirrors + /// Calls the UInt128-based helper; performance mirrors /// the full helper with 5.24–5.95 ns on sparse inputs and 1.99 ns when zeros dominate, but 22.8–59.2 ns on dense blends. /// /// @@ -66,10 +66,10 @@ public class MulMod64Benchmarks /// PrimeSizedModulus 50.490 ns, SparseOperands 5.239 ns, ZeroOperands 1.986 ns. /// [Benchmark] - public ulong GpuCompatibleMulModSimplifiedExtension() + public ulong MulModGpu() { GpuUInt128 gpuUInt128 = new(Input.Left); - return gpuUInt128.MulModSimplified(Input.Right, Input.Modulus); + return gpuUInt128.MulMod(Input.Right, Input.Modulus); } /// diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs new file mode 100644 index 0000000..d4af20d --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs @@ -0,0 +1,94 @@ +using System; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +internal sealed class Pow2MontgomeryModWindowedBenchmarkInputs +{ + public Pow2MontgomeryModWindowedBenchmarkInputs(Pow2MontgomeryModWindowedBenchmarkCase[] cases) + { + Cases = cases; + } + + public Pow2MontgomeryModWindowedBenchmarkCase[] Cases { get; } +} + +internal readonly struct Pow2MontgomeryModWindowedBenchmarkCase +{ + public Pow2MontgomeryModWindowedBenchmarkCase( + ulong exponent, + ulong modulus, + in MontgomeryDivisorData divisor, + ulong cycleLength, + ulong reducedExponent) + { + Exponent = exponent; + Modulus = modulus; + Divisor = divisor; + CycleLength = cycleLength; + ReducedExponent = reducedExponent; + } + + public ulong Exponent { get; } + + public ulong Modulus { get; } + + public MontgomeryDivisorData Divisor { get; } + + public ulong CycleLength { get; } + + public ulong ReducedExponent { get; } +} + +internal static class Pow2MontgomeryModWindowedBenchmarkInputsProvider +{ + private static readonly Lazy Cache = new(Create); + + public static Pow2MontgomeryModWindowedBenchmarkInputs Instance => Cache.Value; + + private static Pow2MontgomeryModWindowedBenchmarkInputs Create() + { + Pow2MontgomeryModWindowedBenchmarkSeed[] seeds = + { + new Pow2MontgomeryModWindowedBenchmarkSeed(138_000_001UL, 1UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(150_000_013UL, 2UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(175_000_019UL, 3UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(190_000_151UL, 4UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(210_000_089UL, 5UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(230_000_039UL, 6UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(260_000_111UL, 7UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(300_000_007UL, 8UL) + }; + + var cases = new Pow2MontgomeryModWindowedBenchmarkCase[seeds.Length]; + for (int i = 0; i < seeds.Length; i++) + { + Pow2MontgomeryModWindowedBenchmarkSeed seed = seeds[i]; + ulong modulus = 2UL * seed.Multiplier * seed.Exponent + 1UL; + MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + ulong cycleLength = MersenneDivisorCycles.CalculateCycleLength(modulus, divisor); + ulong reducedExponent = seed.Exponent % cycleLength; + cases[i] = new Pow2MontgomeryModWindowedBenchmarkCase( + seed.Exponent, + modulus, + divisor, + cycleLength, + reducedExponent); + } + + return new Pow2MontgomeryModWindowedBenchmarkInputs(cases); + } + + private readonly struct Pow2MontgomeryModWindowedBenchmarkSeed + { + public Pow2MontgomeryModWindowedBenchmarkSeed(ulong exponent, ulong multiplier) + { + Exponent = exponent; + Multiplier = multiplier; + } + + public ulong Exponent { get; } + + public ulong Multiplier { get; } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs new file mode 100644 index 0000000..99ad30a --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs @@ -0,0 +1,94 @@ +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +/// +/// Benchmarks CPU Montgomery exponentiation helpers for exponents at or above 138 million. +/// Reuses the shared benchmark input cache so every method processes identical samples. +/// +[MemoryDiagnoser] +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] +public class Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks +{ + private Pow2MontgomeryModWindowedBenchmarkInputs _inputs = null!; + + [GlobalSetup] + public void Setup() + { + _inputs = Pow2MontgomeryModWindowedBenchmarkInputsProvider.Instance; + } + + /// + /// Baseline CPU windowed Montgomery exponentiation that converts the residue back to standard form. + /// + [Benchmark(Baseline = true)] + public ulong WindowedMontgomeryConvert() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedCpu(current.Divisor, keepMontgomery: false); + } + + return checksum; + } + + /// + /// CPU windowed Montgomery exponentiation that keeps the residue in the Montgomery domain. + /// + [Benchmark] + public ulong WindowedMontgomeryKeep() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedCpu(current.Divisor, keepMontgomery: true); + } + + return checksum; + } + + /// + /// Uses the known cycle length to fold the exponent before running the windowed Montgomery ladder. + /// + [Benchmark] + public ulong WindowedMontgomeryWithCycle() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWithCycleCpu(current.CycleLength, current.Divisor); + } + + return checksum; + } + + /// + /// Runs the windowed Montgomery ladder on the reduced exponent that the divisor cycle exposes. + /// + [Benchmark] + public ulong WindowedMontgomeryFromCycleRemainder() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.ReducedExponent.Pow2MontgomeryModFromCycleRemainderCpu(current.Divisor); + } + + return checksum; + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs new file mode 100644 index 0000000..781c33b --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs @@ -0,0 +1,78 @@ +using System; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +/// +/// Benchmarks GPU Montgomery exponentiation helpers for exponents at or above 138 million. +/// The sample set is generated once and reused across all benchmark methods so every variant +/// processes the same data without incurring repeated initialization costs. +/// +[MemoryDiagnoser] +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] +public class Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks +{ + private Pow2MontgomeryModWindowedBenchmarkInputs _inputs = null!; + + [GlobalSetup] + public void Setup() + { + _inputs = Pow2MontgomeryModWindowedBenchmarkInputsProvider.Instance; + } + + /// + /// Baseline GPU windowed Montgomery exponentiation that converts back to the standard residue. + /// + [Benchmark(Baseline = true)] + public ulong WindowedMontgomeryConvert() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedGpu(current.Divisor, keepMontgomery: false); + } + + return checksum; + } + + /// + /// GPU windowed Montgomery exponentiation that keeps the result in the Montgomery domain. + /// + [Benchmark] + public ulong WindowedMontgomeryKeep() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedGpu(current.Divisor, keepMontgomery: true); + } + + return checksum; + } + + /// + /// GPU windowed modular exponentiation without Montgomery reduction for comparison. + /// + [Benchmark] + public ulong WindowedStandardMod() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2ModWindowedGpu(current.Modulus); + } + + return checksum; + } +} diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index 1ef8592..1a8b1ae 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -74,18 +74,6 @@ public class GpuKernelCompilationTests typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetKernel", BindingFlags.NonPublic | BindingFlags.Instance) ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetKernel not found."); - private static readonly MethodInfo DivisorByDivisorGetExponentKernelMethod = - typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetKernelByPrimeExponent", BindingFlags.NonPublic | BindingFlags.Instance) - ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetKernelByPrimeExponent not found."); - - private static readonly MethodInfo DivisorByDivisorGetRemainderDeltaKernelMethod = - typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetRemainderDeltaKernel", BindingFlags.NonPublic | BindingFlags.Instance) - ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetRemainderDeltaKernel not found."); - - private static readonly MethodInfo DivisorByDivisorGetRemainderScanKernelMethod = - typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetRemainderScanKernel", BindingFlags.NonPublic | BindingFlags.Instance) - ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetRemainderScanKernel not found."); - private static readonly MethodInfo DivisorTesterGetKernelMethod = typeof(MersenneNumberDivisorGpuTester).GetMethod("GetKernel", BindingFlags.NonPublic | BindingFlags.Instance) ?? throw new InvalidOperationException("MersenneNumberDivisorGpuTester.GetKernel not found."); @@ -157,9 +145,6 @@ public class GpuKernelCompilationTests yield return Loader("NttGpuMath.ForwardKernel", CompileNttForwardKernel); yield return Loader("NttGpuMath.InverseKernel", CompileNttInverseKernel); yield return Loader("MersenneNumberDivisorByDivisorGpuTester.CheckKernel", CompileDivisorByDivisorCheckKernel); - yield return Loader("MersenneNumberDivisorByDivisorGpuTester.ComputeMontgomeryExponentKernel", CompileDivisorByDivisorExponentKernel); - yield return Loader("MersenneNumberDivisorByDivisorGpuTester.ComputeRemainderDeltasKernel", CompileDivisorByDivisorRemainderDeltaKernel); - yield return Loader("MersenneNumberDivisorByDivisorGpuTester.AccumulateRemaindersKernel", CompileDivisorByDivisorRemainderScanKernel); yield return Loader("MersenneNumberDivisorGpuTester.Kernel", CompileDivisorTesterKernel); yield return Loader("MersenneNumberLucasLehmerGpuTester.Kernel", CompileLucasLehmerKernel); yield return Loader("MersenneNumberLucasLehmerGpuTester.AddSmallKernel", CompileLucasLehmerAddSmallKernel); @@ -281,24 +266,6 @@ public class GpuKernelCompilationTests _ = DivisorByDivisorGetKernelMethod.Invoke(tester, new object[] { accelerator }); } - private static void CompileDivisorByDivisorExponentKernel(Accelerator accelerator) - { - var tester = new MersenneNumberDivisorByDivisorGpuTester(); - _ = DivisorByDivisorGetExponentKernelMethod.Invoke(tester, new object[] { accelerator }); - } - - private static void CompileDivisorByDivisorRemainderDeltaKernel(Accelerator accelerator) - { - var tester = new MersenneNumberDivisorByDivisorGpuTester(); - _ = DivisorByDivisorGetRemainderDeltaKernelMethod.Invoke(tester, new object[] { accelerator }); - } - - private static void CompileDivisorByDivisorRemainderScanKernel(Accelerator accelerator) - { - var tester = new MersenneNumberDivisorByDivisorGpuTester(); - _ = DivisorByDivisorGetRemainderScanKernelMethod.Invoke(tester, new object[] { accelerator }); - } - private static void CompileDivisorTesterKernel(Accelerator accelerator) { var tester = new MersenneNumberDivisorGpuTester(); diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 8f67d06..e19dc8e 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -85,11 +85,17 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_tester_requires_configuration() + public void ByDivisor_tester_defaults_to_zero_limit_before_configuration() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - Action act = () => tester.IsPrime(31UL, out _); - act.Should().Throw(); + + tester.DivisorLimit.Should().Be(0UL); + + ulong[] primes = { 31UL, 37UL }; + ulong[] allowed = new ulong[primes.Length]; + + tester.PrepareCandidates(primes, allowed); + allowed.Should().OnlyContain(value => value == 0UL); } [Fact] @@ -136,6 +142,22 @@ public class MersenneNumberDivisorGpuTesterTests divisorsExhausted.Should().BeTrue(); } + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_session_computes_cycle_when_zero_is_provided() + { + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + tester.ConfigureFromMaxPrime(43UL); + + using var session = tester.CreateDivisorSession(); + ulong[] primes = { 31UL, 37UL, 41UL }; + byte[] hits = new byte[primes.Length]; + + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes, hits); + + hits.Should().ContainInOrder(new byte[] { 0, 1, 0 }); + } + [Fact] [Trait("Category", "Fast")] public void ByDivisor_session_checks_divisors_across_primes() diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs index 10a92a9..2907fe6 100644 --- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs +++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs @@ -42,7 +42,7 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi // Keep these remainder steppers in place so future updates continue reusing the previously computed residues. // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. - var exponentStepper = new ExponentRemainderStepper(cachedData); + var exponentStepper = new ExponentRemainderStepperCpu(cachedData); var cycleStepper = new CycleRemainderStepper(divisorCycle); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 239684b..5ba86d0 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -175,17 +175,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // Keep the divisibility filters aligned with the divisor-cycle generator so the // CPU path never requests cycles that were skipped during cache creation. bool lastIsSeven = (prime & 3UL) == 3UL; + const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); while (divisor.CompareTo(limit) <= 0) { ulong candidate = divisor.Low; processedCount++; - bool admissible = lastIsSeven - ? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9) - : (remainder10 == 1 || remainder10 == 3 || remainder10 == 9); + ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + bool admissible = ((decimalMask >> remainder10) & 1) != 0; - if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0) + if (admissible + && (remainder8 == 1 || remainder8 == 7) + && remainder3 != 0 + && remainder5 != 0 + && remainder7 != 0 + && remainder11 != 0) { MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); ulong divisorCycle; @@ -196,7 +202,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // } // else { - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent( + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( candidate, prime, divisorData, diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index edc92d9..98fe8bf 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -79,7 +79,7 @@ public sealed class DivisorCycleCache throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } - ObserveCycleCacheDivisor(divisor); + // ObserveCycleCacheDivisor(divisor); return cached; } @@ -114,7 +114,7 @@ public sealed class DivisorCycleCache throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } - ObserveCycleCacheDivisor(divisor); + // ObserveCycleCacheDivisor(divisor); cycles[i] = cached; } else diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs similarity index 78% rename from PerfectNumbers.Core/ExponentRemainderStepper.cs rename to PerfectNumbers.Core/ExponentRemainderStepperCpu.cs index 0722880..dc22238 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepper.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs @@ -1,30 +1,28 @@ namespace PerfectNumbers.Core; /// -/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates. +/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path. /// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order. /// -internal struct ExponentRemainderStepper +internal struct ExponentRemainderStepperCpu { private readonly MontgomeryDivisorData _divisor; private readonly ulong _modulus; private readonly ulong _nPrime; private readonly ulong _montgomeryOne; - private ulong _previousExponent; + public ulong PreviousExponent; private ulong _currentMontgomery; - public ExponentRemainderStepper(in MontgomeryDivisorData divisor) + public ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor) { _divisor = divisor; _modulus = divisor.Modulus; _nPrime = divisor.NPrime; _montgomeryOne = divisor.MontgomeryOne; - _previousExponent = 0UL; + PreviousExponent = 0UL; _currentMontgomery = divisor.MontgomeryOne; } - public ulong PreviousExponent => _previousExponent; - public bool MatchesDivisor(in MontgomeryDivisorData divisor) { return _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne; @@ -36,7 +34,7 @@ internal struct ExponentRemainderStepper /// public void Reset() { - _previousExponent = 0UL; + PreviousExponent = 0UL; _currentMontgomery = _montgomeryOne; } @@ -71,38 +69,16 @@ internal struct ExponentRemainderStepper /// The canonical residue produced by 2^exponent mod divisor. public ulong ComputeNextCpu(ulong exponent) { - ulong delta = exponent - _previousExponent; + ulong delta = exponent - PreviousExponent; // TODO: Once divisor cycle lengths are mandatory, pull the delta multiplier from the // single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks. ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); - _previousExponent = exponent; + PreviousExponent = exponent; return ReduceCurrent(); } - /// - /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence. - /// - /// First exponent evaluated on the GPU path. The same divisor metadata must be reused for future calls. - /// The canonical residue produced by 2^exponent mod divisor. - public ulong InitializeGpu(ulong exponent) - { - InitializeGpuState(exponent); - return ReduceCurrent(); - } - - /// - /// Initializes the GPU stepping state and reports whether the seed exponent produces unity. - /// - /// First exponent evaluated on the GPU path. Must match the divisor supplied during construction. - /// when the seed exponent produces a Montgomery residue equal to one. - public bool InitializeGpuIsUnity(ulong exponent) - { - InitializeGpuState(exponent); - return _currentMontgomery == _montgomeryOne; - } - /// /// Advances the CPU stepping state to the next exponent and reports whether it evaluates to unity. /// Callers must seed the stepper through or and then @@ -112,12 +88,12 @@ internal struct ExponentRemainderStepper /// when the canonical residue equals one. public bool ComputeNextIsUnity(ulong exponent) { - ulong delta = exponent - _previousExponent; + ulong delta = exponent - PreviousExponent; // TODO: Reuse the divisor-cycle derived Montgomery delta once the cache exposes single-cycle // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor. ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); - _previousExponent = exponent; + PreviousExponent = exponent; return _currentMontgomery == _montgomeryOne; } @@ -132,7 +108,7 @@ internal struct ExponentRemainderStepper public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity) { _currentMontgomery = montgomeryResult; - _previousExponent = exponent; + PreviousExponent = exponent; isUnity = _currentMontgomery == _montgomeryOne; return true; } @@ -151,7 +127,7 @@ internal struct ExponentRemainderStepper // of relying on the caller-provided delta so incremental scans remain in sync with the snapshot // without computing additional cycles or mutating cache state. _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime); - _previousExponent = exponent; + PreviousExponent = exponent; isUnity = _currentMontgomery == _montgomeryOne; return true; } @@ -159,13 +135,7 @@ internal struct ExponentRemainderStepper private void InitializeCpuState(ulong exponent) { _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _previousExponent = exponent; - } - - private void InitializeGpuState(ulong exponent) - { - _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _previousExponent = exponent; + PreviousExponent = exponent; } private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime); diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs new file mode 100644 index 0000000..d322406 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs @@ -0,0 +1,112 @@ +using System.Runtime.CompilerServices; +using PerfectNumbers.Core; + +namespace PerfectNumbers.Core.Gpu; + +/// +/// Tracks successive residues for a single divisor while kernels advance exponent candidates on the GPU path. +/// Callers must seed the stepper once per divisor and then supply exponents in strictly ascending order. +/// +internal struct ExponentRemainderStepperGpu +{ + private readonly ulong _modulus; + private readonly ReadOnlyGpuUInt128 _modulusWide; + public ulong PreviousExponent; + private GpuUInt128 _currentResidue; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ExponentRemainderStepperGpu(in GpuDivisorPartialData divisor) + { + ulong modulus = divisor.Modulus; + _modulus = modulus; + _modulusWide = divisor.ModulusWide; + PreviousExponent = 0UL; + _currentResidue = GpuUInt128.One; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool MatchesDivisor(in GpuDivisorPartialData divisor) + { + return _modulus == divisor.Modulus; + } + + /// + /// Clears the cached residue so the stepper can be reused for another divisor. + /// Callers must reinitialize the stepper before consuming any residues after a reset. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Reset() + { + PreviousExponent = 0UL; + _currentResidue = GpuUInt128.One; + } + + /// + /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence. + /// The returned value is the canonical residue corresponding to the provided exponent. + /// + /// First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata. + /// The canonical residue produced by 2^exponent mod divisor. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong InitializeGpu(ulong exponent) + { + InitializeStateGpu(exponent); + return _currentResidue.Low; + } + + /// + /// Initializes the GPU stepping state and reports whether the seed exponent already yields unity. + /// + /// First exponent evaluated for the current divisor. Must match the divisor passed to the constructor. + /// when the seed exponent produces a residue equal to one. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool InitializeIsUnityGpu(ulong exponent) + { + InitializeStateGpu(exponent); + return _currentResidue.Low == 1UL; + } + + /// + /// Advances the GPU stepping state to the next exponent in a strictly ascending sequence that shares the same divisor. + /// Callers must invoke beforehand. + /// + /// Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent. + /// The canonical residue produced by 2^exponent mod divisor. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong ComputeNextGpu(ulong exponent) + { + AdvanceStateGpu(exponent); + return _currentResidue.Low; + } + + /// + /// Advances the GPU stepping state to the next exponent and reports whether it evaluates to unity. + /// Callers must seed the stepper through or and then + /// pass exponents in strictly ascending order. + /// + /// Next exponent in the ascending sequence. + /// when the canonical residue equals one. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool ComputeNextIsUnityGpu(ulong exponent) + { + AdvanceStateGpu(exponent); + return _currentResidue.Low == 1UL; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void InitializeStateGpu(ulong exponent) + { + ulong residue = ULongExtensions.Pow2ModWindowedGpuKernel(exponent, _modulus); + _currentResidue = new GpuUInt128(residue); + PreviousExponent = exponent; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void AdvanceStateGpu(ulong exponent) + { + ulong delta = exponent - PreviousExponent; + ulong multiplier = ULongExtensions.Pow2ModWindowedGpuKernel(delta, _modulus); + _currentResidue.MulMod(multiplier, _modulusWide); + PreviousExponent = exponent; + } +} diff --git a/PerfectNumbers.Core/Gpu/GpuDivisorPartialData.cs b/PerfectNumbers.Core/Gpu/GpuDivisorPartialData.cs new file mode 100644 index 0000000..7617ffc --- /dev/null +++ b/PerfectNumbers.Core/Gpu/GpuDivisorPartialData.cs @@ -0,0 +1,23 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core.Gpu; + +internal readonly struct GpuDivisorPartialData +{ + public readonly ulong Modulus; + public readonly ReadOnlyGpuUInt128 ModulusWide; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public GpuDivisorPartialData(ulong modulus) + { + Modulus = modulus; + ModulusWide = new ReadOnlyGpuUInt128(modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public GpuDivisorPartialData(ulong modulus, in ReadOnlyGpuUInt128 modulusWide) + { + Modulus = modulus; + ModulusWide = modulusWide; + } +} diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index af5a386..5159783 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -11,9 +11,7 @@ public struct GpuUInt128 : IComparable, IEquatable private const int NativeModuloChunkBits = 8; private const int NativeModuloChunkBitsMinusOne = NativeModuloChunkBits - 1; - private const int NativeModuloBitMaskTableSize = 1024; private const ulong NativeModuloChunkMask = (1UL << NativeModuloChunkBits) - 1UL; - private static readonly ulong[] NativeModuloBitMasks = CreateNativeModuloBitMasks(); private const int Pow2WindowSizeBits = 8; private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); @@ -902,23 +900,6 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly ulong MulMod(ulong value, ulong modulus) - { - ulong a = Low % modulus; - ulong b = value % modulus; - - if (a == 0UL || b == 0UL) - { - return 0UL; - } - - ulong ulongRange = ulong.MaxValue / a; - return b <= ulongRange - ? (a * b) % modulus - : MulMod64(a, b, modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public readonly ulong MulModSimplified(ulong value, ulong modulus) { ulong modulusLocal = modulus; ulong a = Low % modulusLocal; @@ -1028,10 +1009,10 @@ public struct GpuUInt128 : IComparable, IEquatable private static ulong MultiplyChunkModulo(ulong value, ulong chunk, ulong modulus) { ulong result = 0UL; - ulong[] nativeModuloBitMasks = NativeModuloBitMasks; - for (int bit = 0; bit < NativeModuloChunkBits; bit++) + ulong bitMask = 1UL; + for (int bit = 0; bit < NativeModuloChunkBits; bit++, bitMask <<= 1) { - if ((chunk & nativeModuloBitMasks[bit]) != 0UL) + if ((chunk & bitMask) != 0UL) { result = (result + value) % modulus; } @@ -1548,16 +1529,6 @@ public struct GpuUInt128 : IComparable, IEquatable high = hi; } - private static ulong[] CreateNativeModuloBitMasks() - { - ulong[] masks = new ulong[NativeModuloBitMaskTableSize]; - for (int bit = 0; bit < NativeModuloBitMaskTableSize; bit++) - { - masks[bit] = bit < 64 ? 1UL << bit : 0UL; - } - - return masks; - } } // TODO: Check if the TODO below is still relevant. diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index eece679..accaa5f 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -7,280 +7,82 @@ namespace PerfectNumbers.Core.Gpu; internal static class DivisorByDivisorKernels { - public static void CheckKernel(Index1D index, ArrayView divisors, ArrayView exponents, ArrayView hits) - { - MontgomeryDivisorData divisor = divisors[index]; - ulong modulus = divisor.Modulus; - // The GPU by-divisor pipeline only materializes odd moduli greater than one (q = 2kp + 1), - // so the defensive guard for invalid values stays disabled here to keep the kernel branch-free. - ulong exponent = exponents[index]; - ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent); - hits[index] = montgomeryResult == 1UL ? (byte)1 : (byte)0; - } - - public static void GenerateCandidatesAndGapsKernel( - Index1D index, - ArrayView candidates, - ArrayView gaps, - ArrayView nextStarts, - GpuUInt128 startValue, - GpuUInt128 stride) - { - int globalIndex = index; - // EvenPerfectBitScanner always launches this kernel with an extent matching the candidate view length, - // so the bounds guard stays commented out to keep the hot path branch-free. - // int candidateLength = (int)candidates.Length; - // if (globalIndex >= candidateLength) - // { - // return; - // } - - ulong offset = (ulong)globalIndex; - GpuUInt128 candidateValue = stride; - candidateValue.Mul(offset); - candidateValue.Add(startValue); - candidates[globalIndex] = candidateValue.Low; - - // The scanner supplies gap buffers that mirror the candidate view length, - // so the additional bounds guard remains commented out for the same reason as above. - // int gapLength = (int)gaps.Length; - // if (globalIndex >= gapLength) - // { - // return; - // } - - ulong previousOffset = offset > 0UL ? offset - 1UL : 0UL; - GpuUInt128 previousValue = stride; - previousValue.Mul(previousOffset); - previousValue.Add(startValue); - - GpuUInt128 gapValue = candidateValue; - gapValue.Sub(previousValue); - gaps[globalIndex] = gapValue.Low; - - if (globalIndex == 0) - { - // EvenPerfectBitScanner provides a single-element nextStarts view per launch, so the guard stays commented out. - // int nextStartLength = (int)nextStarts.Length; - // if (nextStartLength <= 0) - // { - // return; - // } - - ulong candidateCount = (ulong)candidates.Length; - GpuUInt128 nextValue = stride; - nextValue.Mul(candidateCount); - nextValue.Add(startValue); - nextStarts[0] = nextValue; - } - } - - public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) - { - deltas[index] = (byte)(gaps[index] % modulus); - } - - public static void AccumulateRemaindersKernel( - Index1D index, - ArrayView deltas, - ArrayView remainders, - byte baseRemainder, - byte modulus) - { - int globalIndex = index; - // EvenPerfectBitScanner launches this stream kernel with an extent matching the remainder span length, - // so the defensive guard stays commented out to keep the kernel branch-free. - // int length = (int)remainders.Length; - // if (globalIndex >= length) - // { - // return; - // } - - int localIndex = Group.IdxX; - int groupSize = Group.Dimension.X; - var shared = SharedMemory.GetDynamic(); - - shared[localIndex] = deltas[globalIndex]; - Group.Barrier(); - - int offset = 1; - int mod = modulus; - while (offset < groupSize) - { - bool apply = localIndex >= offset; - int addend = apply ? shared[localIndex - offset] : 0; - - Group.Barrier(); - - int current = shared[localIndex]; - int sum = apply ? current + addend : current; - int wrap = sum >= mod ? 1 : 0; - sum -= wrap * mod; - shared[localIndex] = apply ? sum : current; - - Group.Barrier(); - - offset <<= 1; - } - - int remainder = baseRemainder + shared[localIndex]; - int remainderWrap = remainder >= mod ? 1 : 0; - remainder -= remainderWrap * mod; - remainders[globalIndex] = (byte)remainder; - } - - public static void EvaluateCandidateMaskKernel( - Index1D index, - ArrayView remainder10, - ArrayView remainder8, - ArrayView remainder5, - ArrayView remainder3, - byte lastIsSevenFlag, - ArrayView mask) - { - int globalIndex = index; - // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the mask length, - // so the defensive guard stays commented out to keep the kernel branch-free. - // int length = (int)mask.Length; - // if (globalIndex >= length) - // { - // return; - // } - - byte value10 = remainder10[globalIndex]; - bool lastIsSeven = lastIsSevenFlag != 0; - bool acceptSeven = value10 == 3 || value10 == 7 || value10 == 9; - bool acceptNonSeven = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9; - bool accept10 = lastIsSeven ? acceptSeven : acceptNonSeven; - byte value8 = remainder8[globalIndex]; - bool accept8 = value8 == 1 || value8 == 7; - bool accept5 = remainder5[globalIndex] != 0; - bool accept3 = remainder3[globalIndex] != 0; - byte accepted = (byte)((accept10 && accept8 && accept3 && accept5) ? 1 : 0); - mask[globalIndex] = accepted; - } - - public static void CompactCandidateMaskKernel( - Index1D index, - ArrayView mask, - ArrayView sourceCandidates, - ArrayView compactedCandidates, - ArrayView compactedCount) - { - int globalIndex = index; - int length = (int)mask.Length; - // EvenPerfectBitScanner launches this stream kernel with an extent matching the mask length, - // so the defensive guard stays commented out to keep the kernel branch-free. - // if (globalIndex >= length) - // { - // return; - // } - - int localIndex = Group.IdxX; - int groupSize = Group.Dimension.X; - var shared = SharedMemory.GetDynamic(); - - int prefixBase = 0; - int mappingBase = length; - int fallbackBase = length * 2; - // The dynamic shared-memory slice stores prefix sums, compacted index mappings, and fallback slots - // in separate regions so the scatter stage can avoid divergent branches. - - byte maskValue = mask[globalIndex]; - int accepted = maskValue != 0 ? 1 : 0; - shared[prefixBase + localIndex] = accepted; - shared[fallbackBase + localIndex] = globalIndex; - - Group.Barrier(); - - int offset = 1; - while (offset < groupSize) - { - int addendIndex = prefixBase + localIndex - offset; - int addend = localIndex >= offset ? shared[addendIndex] : 0; - - Group.Barrier(); - - int currentIndex = prefixBase + localIndex; - int current = shared[currentIndex]; - int sum = current + addend; - int shouldUpdate = localIndex >= offset ? 1 : 0; - shared[currentIndex] = shouldUpdate != 0 ? sum : current; - - Group.Barrier(); - - offset <<= 1; - } - - int inclusivePrefix = shared[prefixBase + localIndex]; - int exclusiveIndex = inclusivePrefix - accepted; - int targetBase = accepted != 0 ? mappingBase : fallbackBase; - int targetOffset = accepted != 0 ? exclusiveIndex : localIndex; - shared[targetBase + targetOffset] = globalIndex; - - Group.Barrier(); - - int totalAccepted = shared[prefixBase + length - 1]; - int writerMask = localIndex < totalAccepted ? 1 : 0; - int lookupBase = writerMask != 0 ? mappingBase : fallbackBase; - int sourceIndex = shared[lookupBase + localIndex]; - - ulong candidateValue = sourceCandidates[sourceIndex]; - int destinationIndex = localIndex; - ulong existingValue = compactedCandidates[destinationIndex]; - compactedCandidates[destinationIndex] = writerMask != 0 ? candidateValue : existingValue; - - int lastLaneMask = globalIndex == length - 1 ? 1 : 0; - int countValue = shared[prefixBase + localIndex]; - if (lastLaneMask != 0) - { - compactedCount[0] = countValue; - } - } - - public static void PublishFirstHitKernel( + public static void CheckKernel( Index1D index, + ArrayView divisors, + ArrayView exponentOffsets, + ArrayView exponentCounts, + ArrayView exponents, + ArrayView divisorCycles, + ArrayView firstCycleRemainders, ArrayView hits, ArrayView firstHit) { int globalIndex = index; - // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the hit view length, - // so the defensive guard stays commented out to keep the hot path branch-free. - // int length = (int)hits.Length; - // if (globalIndex >= length) - // { - // return; - // } + int count = exponentCounts[globalIndex]; + if (count <= 0) + { + return; + } - byte hitValue = hits[globalIndex]; - int candidate = hitValue != 0 ? globalIndex : int.MaxValue; - Atomic.Min(ref firstHit[0], candidate); - } + int offset = exponentOffsets[globalIndex]; + GpuDivisorPartialData divisor = divisors[globalIndex]; + ulong cycleLength = divisorCycles[globalIndex]; + ulong firstCycleRemainder = firstCycleRemainders[globalIndex]; - public static void ConvertMontgomeryResultsToHitsKernel( - Index1D index, - ArrayView results, - ArrayView hits) - { - int globalIndex = index; - // EvenPerfectBitScanner launches this kernel with an extent matching the result and hit spans, - // so the defensive bounds guard stays commented out to keep the conversion branch-free. - // int length = (int)results.Length; - // if (globalIndex >= length) - // { - // return; - // } + var stepper = new ExponentRemainderStepperGpu(divisor); - ulong montgomeryResult = results[globalIndex]; - hits[globalIndex] = montgomeryResult == 1UL ? (byte)1 : (byte)0; - } + ulong firstExponent = exponents[offset]; + bool firstUnity = stepper.InitializeIsUnityGpu(firstExponent); - public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) - { - ulong modulus = divisor.Modulus; - // Each divisor flowing through this kernel follows q = 2kp + 1 with k >= 1, so modulus is always odd and exceeds one. - ulong exponent = exponents[index]; - ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponent); - results[index] = montgomeryResult; + byte firstHitValue; + ulong previousExponent = firstExponent; + ulong cycleRemainder = 0UL; + if (cycleLength == 0UL) + { + firstHitValue = firstUnity ? (byte)1 : (byte)0; + } + else + { + cycleRemainder = firstCycleRemainder; + firstHitValue = cycleRemainder == 0UL && firstUnity ? (byte)1 : (byte)0; + } + + hits[offset] = firstHitValue; + byte anyHit = firstHitValue; + + for (int i = 1; i < count; i++) + { + ulong exponent = exponents[offset + i]; + ulong delta = exponent - previousExponent; + previousExponent = exponent; + + if (cycleLength != 0UL) + { + var remainderValue = new GpuUInt128(cycleRemainder); + remainderValue.AddMod(delta, cycleLength); + cycleRemainder = remainderValue.Low; + + byte hit = 0; + if (cycleRemainder == 0UL) + { + hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; + } + + hits[offset + i] = hit; + anyHit |= hit; + } + else + { + byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; + hits[offset + i] = hit; + anyHit |= hit; + } + } + + if (anyHit != 0 && firstHit.Length > 0) + { + Atomic.Min(ref firstHit[0], globalIndex); + } } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 33ed320..1a2e728 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,3 +1,4 @@ +using System; using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; @@ -10,12 +11,11 @@ using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu; /// -/// Scans Mersenne divisors on the GPU for prime exponents p >= 31 using per-divisor Montgomery data. +/// Scans Mersenne divisors on the GPU for prime exponents p >= 31 using cached GPU divisor partial data. +/// Consumers must call before invoking other members so the divisor limits are populated. /// public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDivisorByDivisorTester { - private const int RemainderSlotCount = 4; - private int _gpuBatchSize = GpuConstants.ScanBatchSize; // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, // so the synchronization fields from the previous implementation remain commented out here. @@ -23,57 +23,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private ulong _divisorLimit; // private bool _isConfigured; - private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); - private readonly ConcurrentDictionary, ArrayView, byte>> _remainderDeltaKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, byte, byte>> _remainderScanKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView, byte, ArrayView>> _candidateMaskKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView>> _candidateCompactionKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView>> _hitDetectionKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView>> _hitConversionKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>> _kernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); - private Action, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(DivisorByDivisorKernels.CheckKernel)); + private Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => + acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(DivisorByDivisorKernels.CheckKernel)); - private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.ComputeMontgomeryExponentKernel)); - - private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) => - _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(DivisorByDivisorKernels.ComputeRemainderDeltasKernel)); - - private Action, ArrayView, byte, byte> GetRemainderScanKernel(Accelerator accelerator) => - _remainderScanKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadStreamKernel, ArrayView, byte, byte>(DivisorByDivisorKernels.AccumulateRemaindersKernel)); - - private Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> GetCandidateMaskKernel(Accelerator accelerator) => - _candidateMaskKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); - - private Action, ArrayView, ArrayView, ArrayView> GetCandidateCompactionKernel(Accelerator accelerator) => - _candidateCompactionKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadStreamKernel, ArrayView, ArrayView, ArrayView>(DivisorByDivisorKernels.CompactCandidateMaskKernel)); - - private Action, ArrayView> GetHitDetectionKernel(Accelerator accelerator) => - _hitDetectionKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.PublishFirstHitKernel)); - - private Action, ArrayView> GetResultConversionKernel(Accelerator accelerator) => - _hitConversionKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.ConvertMontgomeryResultsToHitsKernel)); - - private Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => - _candidateGenerationKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, GpuUInt128, GpuUInt128>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel)); public int GpuBatchSize { @@ -144,30 +102,22 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv composite = CheckDivisors( prime, allowedMax, - accelerator, kernel, - resources.DivisorsBuffer, + resources.DivisorDataBuffer, + resources.OffsetBuffer, + resources.CountBuffer, + resources.CycleBuffer, + resources.FirstRemainderBuffer, resources.ExponentBuffer, resources.HitsBuffer, resources.HitIndexBuffer, - resources.CandidateBuffer, - resources.NextStartBuffer, - resources.DivisorDeltaBuffer, - resources.RemainderDeltaBuffer, - resources.RemainderBuffer, - resources.FilterCountBuffer, - resources.RemainderDeltaKernel, - resources.RemainderScanKernel, - resources.CandidateMaskKernel, - resources.CandidateCompactionKernel, - resources.CandidateGenerationKernel, - resources.HitDetectionKernel, resources.Divisors, resources.Exponents, resources.DivisorData, - resources.CycleCandidates, - resources.CycleLengths, - resources.CycleCapacity, + resources.Offsets, + resources.Counts, + resources.Cycles, + resources.FirstRemainders, out lastProcessed, out coveredRange, out processedCount); @@ -218,46 +168,33 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static bool CheckDivisors( ulong prime, ulong allowedMax, - Accelerator accelerator, - Action, ArrayView, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, + Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> kernel, + MemoryBuffer1D divisorDataBuffer, + MemoryBuffer1D offsetBuffer, + MemoryBuffer1D countBuffer, + MemoryBuffer1D cycleBuffer, + MemoryBuffer1D firstRemainderBuffer, MemoryBuffer1D exponentBuffer, MemoryBuffer1D hitsBuffer, MemoryBuffer1D hitIndexBuffer, - MemoryBuffer1D candidateBuffer, - MemoryBuffer1D nextStartBuffer, - MemoryBuffer1D divisorDeltaBuffer, - MemoryBuffer1D remainderDeltaBuffer, - MemoryBuffer1D remainderBuffer, - MemoryBuffer1D filterCountBuffer, - Action, ArrayView, byte> remainderDeltaKernel, - Action, ArrayView, byte, byte> remainderScanKernel, - Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, - Action, ArrayView, ArrayView, ArrayView> candidateCompactionKernel, - Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, - Action, ArrayView> hitDetectionKernel, in ulong[] divisors, ulong[] exponents, - MontgomeryDivisorData[] divisorData, - ulong[] cycleCandidates, - ulong[] cycleLengths, - int cycleCapacity, + GpuDivisorPartialData[] divisorData, + int[] offsets, + int[] counts, + ulong[] cycles, + ulong[] firstRemainders, out ulong lastProcessed, out bool coveredRange, out ulong processedCount) { - int batchCapacity = (int)divisorsBuffer.Length; + int batchCapacity = (int)divisorDataBuffer.Length; bool composite = false; bool processedAll = false; processedCount = 0UL; lastProcessed = 0UL; - cycleCapacity = Math.Max(1, cycleCapacity); - int chunkCapacity = Math.Min(batchCapacity, cycleCapacity); - int maxThreadsPerGroup = (int)accelerator.MaxNumThreadsPerGroup; - bool clampThreads = maxThreadsPerGroup > 0; - chunkCapacity = clampThreads && chunkCapacity > maxThreadsPerGroup ? maxThreadsPerGroup : chunkCapacity; - chunkCapacity = chunkCapacity <= 0 ? cycleCapacity : chunkCapacity; + int chunkCapacity = Math.Max(1, batchCapacity); UInt128 twoP128 = (UInt128)prime << 1; UInt128 allowedMax128 = allowedMax; @@ -280,111 +217,86 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv byte step8 = (byte)(twoP128 % 8UL); byte step5 = (byte)(twoP128 % 5UL); byte step3 = (byte)(twoP128 % 3UL); + byte step7 = (byte)(twoP128 % 7UL); + byte step11 = (byte)(twoP128 % 11UL); UInt128 currentDivisor128 = firstDivisor128; byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL)); byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL)); byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL)); byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL)); + byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL)); + byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL)); bool lastIsSeven = (prime & 3UL) == 3UL; - byte lastIsSevenFlag = lastIsSeven ? (byte)1 : (byte)0; - DivisorCycleCache cycleCache = DivisorCycleCache.Shared; - int remainderStride = (int)divisorDeltaBuffer.Length; - remainderStride = remainderStride <= 0 ? chunkCapacity : remainderStride; ArrayPool ulongPool = ThreadStaticPools.UlongPool; ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity); - var viewCache = new BufferSliceCache(divisorDeltaBuffer, remainderDeltaBuffer, remainderBuffer, hitsBuffer, remainderStride); - Span candidateStorage = cycleCandidates.AsSpan(); Span filteredStorage = filteredDivisors.AsSpan(); - Span cycleLengthStorage = cycleLengths.AsSpan(); - Span divisorDataStorage = divisorData.AsSpan(); Span divisorStorage = divisors.AsSpan(); Span exponentStorage = exponents.AsSpan(); + Span divisorDataStorage = divisorData.AsSpan(); + Span offsetStorage = offsets.AsSpan(); + Span countStorage = counts.AsSpan(); + Span cycleStorage = cycles.AsSpan(); + Span remainderStorage = firstRemainders.AsSpan(); while (currentK <= maxK && !composite) { int chunkCount = Math.Min(chunkCapacity, batchCapacity); ulong remainingK = maxK - currentK + 1UL; - chunkCount = (ulong)chunkCount > remainingK ? (int)remainingK : chunkCount; - bool hasChunk = chunkCount > 0; - int effectiveCount = hasChunk ? chunkCount : 0; - if (!hasChunk) + if ((ulong)chunkCount > remainingK) + { + chunkCount = (int)remainingK; + } + + if (chunkCount <= 0) { processedAll = true; break; } - Span candidateSpan = candidateStorage[..effectiveCount]; + int filteredCount = 0; + UInt128 nextDivisor128 = currentDivisor128; - viewCache.Ensure(effectiveCount); + byte localRemainder10 = remainder10; + byte localRemainder8 = remainder8; + byte localRemainder5 = remainder5; + byte localRemainder3 = remainder3; + byte localRemainder7 = remainder7; + byte localRemainder11 = remainder11; - ArrayView1D gapView = viewCache.DivisorDelta; - ArrayView1D deltaView = viewCache.RemainderDelta; - ArrayView1D remainder10View = viewCache.Remainder10; - ArrayView1D remainder8View = viewCache.Remainder8; - ArrayView1D remainder5View = viewCache.Remainder5; - ArrayView1D remainder3View = viewCache.Remainder3; - ArrayView1D maskView = viewCache.Hits; - ArrayView1D filteredCountView = filterCountBuffer.View; - - UInt128 localDivisor = currentDivisor128; - - ArrayView1D candidateView = candidateBuffer.View.SubView(0, effectiveCount); - ArrayView1D nextStartView = nextStartBuffer.View.SubView(0, 1); - GpuUInt128 gpuStartValue = new GpuUInt128(localDivisor); - GpuUInt128 gpuStride = new GpuUInt128(twoP128); - GpuUInt128 nextDivisorGpu = GenerateCandidatesOnGpu( - effectiveCount, - gpuStartValue, - gpuStride, - candidateView, - gapView, - nextStartView, - candidateGenerationKernel); - if (effectiveCount > 0) + for (int i = 0; i < chunkCount; i++) { - candidateView.CopyToCPU(ref MemoryMarshal.GetReference(candidateSpan), effectiveCount); - processedCount += (ulong)effectiveCount; - lastProcessed = candidateSpan[effectiveCount - 1]; - currentDivisor128 = (UInt128)nextDivisorGpu; + ulong candidate = (ulong)nextDivisor128; + + if (CandidatePassesHeuristics(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastIsSeven)) + { + filteredStorage[filteredCount++] = candidate; + } + + nextDivisor128 += twoP128; + localRemainder10 = AddMod(localRemainder10, step10, 10); + localRemainder8 = AddMod(localRemainder8, step8, 8); + localRemainder5 = AddMod(localRemainder5, step5, 5); + localRemainder3 = AddMod(localRemainder3, step3, 3); + localRemainder7 = AddMod(localRemainder7, step7, 7); + localRemainder11 = AddMod(localRemainder11, step11, 11); } - ComputeRemaindersOnGpu(effectiveCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(effectiveCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(effectiveCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(effectiveCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel); + remainder10 = localRemainder10; + remainder8 = localRemainder8; + remainder5 = localRemainder5; + remainder3 = localRemainder3; + remainder7 = localRemainder7; + remainder11 = localRemainder11; - BuildCandidateMaskOnGpu(effectiveCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, candidateMaskKernel); + processedCount += (ulong)chunkCount; - byte lastRemainder10 = effectiveCount > 0 ? FetchLastRemainder(remainder10View, effectiveCount) : (byte)0; - byte lastRemainder8 = effectiveCount > 0 ? FetchLastRemainder(remainder8View, effectiveCount) : (byte)0; - byte lastRemainder5 = effectiveCount > 0 ? FetchLastRemainder(remainder5View, effectiveCount) : (byte)0; - byte lastRemainder3 = effectiveCount > 0 ? FetchLastRemainder(remainder3View, effectiveCount) : (byte)0; + UInt128 lastDivisor128 = nextDivisor128 - twoP128; + lastProcessed = (ulong)lastDivisor128; + currentDivisor128 = nextDivisor128; - remainder10 = AddModGpu(lastRemainder10, step10, 10); - remainder8 = AddModGpu(lastRemainder8, step8, 8); - remainder5 = AddModGpu(lastRemainder5, step5, 5); - remainder3 = AddModGpu(lastRemainder3, step3, 3); - - int filteredCount = CompactCandidatesOnGpu( - effectiveCount, - maskView, - candidateView, - candidateView, - filteredCountView, - candidateCompactionKernel); - - if (filteredCount > 0) - { - Span filteredSlice = filteredStorage[..filteredCount]; - ArrayView1D filteredView = candidateView.SubView(0, filteredCount); - // The compaction kernel rewrites the accepted divisors into the front of the candidate buffer, - // so copying the compacted prefix is enough to build the CPU span for cycle checks. - filteredView.CopyToCPU(ref MemoryMarshal.GetReference(filteredSlice), filteredCount); - } - - currentK += (ulong)effectiveCount; + currentK += (ulong)chunkCount; processedAll = currentK > maxK; if (filteredCount == 0) @@ -393,15 +305,37 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } Span filteredDivisorsSpan = filteredStorage[..filteredCount]; - Span cycleSpan = cycleLengthStorage[..filteredCount]; - cycleCache.GetCycleLengths(filteredDivisorsSpan, cycleSpan); + + Span divisorSpan = divisorStorage; + Span exponentSpan = exponentStorage; + Span divisorDataSpan = divisorDataStorage; + Span offsetSpan = offsetStorage; + Span countSpan = countStorage; + Span cycleSpan = cycleStorage; + Span remainderSpan = remainderStorage; int admissibleCount = 0; + int exponentIndex = 0; for (int i = 0; i < filteredCount; i++) { - int add = cycleSpan[i] == prime ? 1 : 0; - filteredDivisorsSpan[admissibleCount] = filteredDivisorsSpan[i]; - admissibleCount += add; + ulong divisorValue = filteredDivisorsSpan[i]; + MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); + ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData); + if (divisorCycle != prime) + { + continue; + } + + divisorSpan[admissibleCount] = divisorValue; + exponentSpan[exponentIndex] = prime; + divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); + offsetSpan[admissibleCount] = exponentIndex; + countSpan[admissibleCount] = 1; + cycleSpan[admissibleCount] = divisorCycle; + remainderSpan[admissibleCount] = 0UL; + + admissibleCount++; + exponentIndex++; } if (admissibleCount == 0) @@ -409,27 +343,30 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv continue; } - Span divisorDataSpan = divisorDataStorage[..admissibleCount]; - Span divisorSpan = divisorStorage[..admissibleCount]; - Span exponentSpan = exponentStorage[..admissibleCount]; - ArrayView1D divisorView = divisorsBuffer.View.SubView(0, admissibleCount); - ArrayView1D exponentView = exponentBuffer.View.SubView(0, admissibleCount); - ArrayView1D hitsView = hitsBuffer.View.SubView(0, admissibleCount); + ArrayView1D divisorView = divisorDataBuffer.View.SubView(0, admissibleCount); + ArrayView1D offsetView = offsetBuffer.View.SubView(0, admissibleCount); + ArrayView1D countView = countBuffer.View.SubView(0, admissibleCount); + ArrayView1D cycleView = cycleBuffer.View.SubView(0, admissibleCount); + ArrayView1D remainderView = firstRemainderBuffer.View.SubView(0, admissibleCount); + ArrayView1D exponentView = exponentBuffer.View.SubView(0, exponentIndex); + ArrayView1D hitsView = hitsBuffer.View.SubView(0, exponentIndex); ArrayView1D hitIndexView = hitIndexBuffer.View.SubView(0, 1); - for (int i = 0; i < admissibleCount; i++) - { - ulong divisorValue = filteredDivisorsSpan[i]; - divisorSpan[i] = divisorValue; - exponentSpan[i] = prime; - divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorValue); - } - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), admissibleCount); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), admissibleCount); - kernel(admissibleCount, divisorView, exponentView, hitsView); + offsetView.CopyFromCPU(ref MemoryMarshal.GetReference(offsetSpan), admissibleCount); + countView.CopyFromCPU(ref MemoryMarshal.GetReference(countSpan), admissibleCount); + cycleView.CopyFromCPU(ref MemoryMarshal.GetReference(cycleSpan), admissibleCount); + remainderView.CopyFromCPU(ref MemoryMarshal.GetReference(remainderSpan), admissibleCount); + exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), exponentIndex); - int hitIndex = FindFirstHitOnGpu(admissibleCount, hitsView, hitIndexView, hitDetectionKernel); + int sentinel = int.MaxValue; + hitIndexView.CopyFromCPU(ref sentinel, 1); + + kernel(admissibleCount, divisorView, offsetView, countView, exponentView, cycleView, remainderView, hitsView, hitIndexView); + + int firstHit = sentinel; + hitIndexView.CopyToCPU(ref firstHit, 1); + int hitIndex = firstHit >= admissibleCount ? -1 : firstHit; bool hitFound = hitIndex >= 0; composite = hitFound; int lastIndex = admissibleCount - 1; @@ -442,218 +379,70 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return composite; } - private ref struct BufferSliceCache - { - private readonly MemoryBuffer1D _divisorDeltaBuffer; - private readonly MemoryBuffer1D _remainderDeltaBuffer; - private readonly MemoryBuffer1D _remainderBuffer; - private readonly MemoryBuffer1D _hitsBuffer; - private readonly int _remainderStride; - private int _length; - private ArrayView1D _divisorDelta; - private ArrayView1D _remainderDelta; - private ArrayView1D _remainder10; - private ArrayView1D _remainder8; - private ArrayView1D _remainder5; - private ArrayView1D _remainder3; - private ArrayView1D _hits; - internal BufferSliceCache( - MemoryBuffer1D divisorDeltaBuffer, - MemoryBuffer1D remainderDeltaBuffer, - MemoryBuffer1D remainderBuffer, - MemoryBuffer1D hitsBuffer, - int remainderStride) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ResolveDivisorCycle(ulong divisor, ulong prime, in MontgomeryDivisorData divisorData) + { + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL) { - _divisorDeltaBuffer = divisorDeltaBuffer; - _remainderDeltaBuffer = remainderDeltaBuffer; - _remainderBuffer = remainderBuffer; - _hitsBuffer = hitsBuffer; - _remainderStride = remainderStride; - _length = -1; - _divisorDelta = default; - _remainderDelta = default; - _remainder10 = default; - _remainder8 = default; - _remainder5 = default; - _remainder3 = default; - _hits = default; + return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); } - internal void Ensure(int length) - { - if (_length == length) - { - return; - } + return computedCycle; + } - _length = length; - ArrayView1D remainderBaseView = _remainderBuffer.View; - _divisorDelta = _divisorDeltaBuffer.View.SubView(0, length); - _remainderDelta = _remainderDeltaBuffer.View.SubView(0, length); - _remainder10 = remainderBaseView.SubView(0, length); - _remainder8 = remainderBaseView.SubView(_remainderStride, length); - _remainder5 = remainderBaseView.SubView(_remainderStride * 2, length); - _remainder3 = remainderBaseView.SubView(_remainderStride * 3, length); - _hits = _hitsBuffer.View.SubView(0, length); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool CandidatePassesHeuristics( + byte remainder10, + byte remainder8, + byte remainder3, + byte remainder5, + byte remainder7, + byte remainder11, + bool lastIsSeven) + { + const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + + ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + if (((decimalMask >> remainder10) & 1) == 0) + { + return false; } - internal ArrayView1D DivisorDelta => _divisorDelta; + if (remainder8 != 1 && remainder8 != 7) + { + return false; + } - internal ArrayView1D RemainderDelta => _remainderDelta; + if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0) + { + return false; + } - internal ArrayView1D Remainder10 => _remainder10; - - internal ArrayView1D Remainder8 => _remainder8; - - internal ArrayView1D Remainder5 => _remainder5; - - internal ArrayView1D Remainder3 => _remainder3; - - internal ArrayView1D Hits => _hits; + return true; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static GpuUInt128 GenerateCandidatesOnGpu( - int count, - GpuUInt128 startValue, - GpuUInt128 stride, - ArrayView1D candidateView, - ArrayView1D gapView, - ArrayView1D nextStartView, - Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> kernel) - { - // EvenPerfectBitScanner always launches this helper with a positive count, - // so the guard stays commented out to keep the hot path branch-free. - // if (count <= 0) - // { - // return GpuUInt128.Zero; - // } - - kernel(count, candidateView, gapView, nextStartView, startValue, stride); - - GpuUInt128 nextStart = default; - nextStartView.SubView(0, 1).CopyToCPU(ref nextStart, 1); - return nextStart; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void ComputeRemaindersOnGpu( - int count, - byte baseRemainder, - byte modulus, - ArrayView1D gapView, - ArrayView1D deltaView, - ArrayView1D remainderView, - Action, ArrayView, byte> deltaKernel, - Action, ArrayView, byte, byte> scanKernel) - { - // EvenPerfectBitScanner always supplies a positive count here, - // so the guard stays commented out to keep the hot path branch-free. - // if (count <= 0) - // { - // return; - // } - - deltaKernel(count, gapView, deltaView, modulus); - SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); - Index1D gridDim = new Index1D(1); - Index1D groupDim = new Index1D(count); - KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig); - - scanKernel(kernelConfig, new Index1D(count), deltaView, remainderView, baseRemainder, modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void BuildCandidateMaskOnGpu( - int count, - ArrayView1D remainder10View, - ArrayView1D remainder8View, - ArrayView1D remainder5View, - ArrayView1D remainder3View, - byte lastIsSevenFlag, - ArrayView1D maskView, - Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> maskKernel) - { - // EvenPerfectBitScanner always supplies a positive count here, - // so the guard stays commented out to keep the hot path branch-free. - // if (count <= 0) - // { - // return; - // } - - maskKernel(count, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView); - } - - private static int CompactCandidatesOnGpu( - int count, - ArrayView1D maskView, - ArrayView1D sourceView, - ArrayView1D destinationView, - ArrayView1D countView, - Action, ArrayView, ArrayView, ArrayView> compactionKernel) - { - // EvenPerfectBitScanner always compacts with a positive count, - // so the guard stays commented out to keep the hot path branch-free. - // if (count <= 0) - // { - // return 0; - // } - - SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count * 3); - Index1D gridDim = new Index1D(1); - Index1D groupDim = new Index1D(count); - KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig); - - compactionKernel(kernelConfig, new Index1D(count), maskView, sourceView, destinationView, countView); - - int filteredCount = 0; - countView.SubView(0, 1).CopyToCPU(ref filteredCount, 1); - return filteredCount; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int FindFirstHitOnGpu( - int count, - ArrayView1D hitsView, - ArrayView1D hitIndexView, - Action, ArrayView> hitKernel) - { - // EvenPerfectBitScanner only calls this helper with positive counts, - // so the guard stays commented out to keep the hot path branch-free. - // if (count <= 0) - // { - // return -1; - // } - - int sentinel = int.MaxValue; - hitIndexView.SubView(0, 1).CopyFromCPU(ref sentinel, 1); - - hitKernel(count, hitsView, hitIndexView); - - int firstHit = sentinel; - hitIndexView.SubView(0, 1).CopyToCPU(ref firstHit, 1); - return firstHit >= count ? -1 : firstHit; - } - - private static byte FetchLastRemainder(ArrayView1D remainderView, int count) - { - // Callers guard against zero counts before invoking this helper to keep the hot path branch-free. - byte value = 0; - remainderView.SubView(count - 1, 1).CopyToCPU(ref value, 1); - return value; - } - - private static byte AddModGpu(byte value, byte delta, byte modulus) + private static byte AddMod(byte value, byte delta, byte modulus) { int result = value + delta; - int firstWrap = result >= modulus ? 1 : 0; - result -= firstWrap * modulus; - int secondWrap = result >= modulus ? 1 : 0; - result -= secondWrap * modulus; + if (result >= modulus) + { + result -= modulus; + if (result >= modulus) + { + result -= modulus; + } + } + return (byte)result; } + + + + private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) { // The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1, @@ -691,10 +480,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; - private Action, ArrayView> _kernel = null!; - private Action, ArrayView> _hitConversionKernel = null!; + private Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> _kernel = null!; + private MemoryBuffer1D _divisorBuffer = null!; + private MemoryBuffer1D _offsetBuffer = null!; + private MemoryBuffer1D _countBuffer = null!; + private MemoryBuffer1D _cycleBuffer = null!; + private MemoryBuffer1D _firstRemainderBuffer = null!; + private MemoryBuffer1D _firstHitBuffer = null!; private MemoryBuffer1D _exponentsBuffer = null!; - private MemoryBuffer1D _resultsBuffer = null!; private MemoryBuffer1D _hitBuffer = null!; private ulong[] _hostBuffer = null!; private int _capacity; @@ -713,14 +506,47 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _disposed = false; } - private void EnsureExecutionResourcesLocked() + private void EnsureExecutionResourcesLocked(int requiredCapacity) { - _kernel = _owner.GetKernelByPrimeExponent(_accelerator); - _hitConversionKernel = _owner.GetResultConversionKernel(_accelerator); - _exponentsBuffer = _accelerator.Allocate1D(_capacity); - _resultsBuffer = _accelerator.Allocate1D(_capacity); - _hitBuffer = _accelerator.Allocate1D(_capacity); - _hostBuffer = ThreadStaticPools.UlongPool.Rent(_capacity); + if (_kernel == null) + { + _kernel = _owner.GetKernel(_accelerator); + } + + if (_divisorBuffer == null) + { + _divisorBuffer = _accelerator.Allocate1D(1); + _offsetBuffer = _accelerator.Allocate1D(1); + _countBuffer = _accelerator.Allocate1D(1); + _cycleBuffer = _accelerator.Allocate1D(1); + _firstRemainderBuffer = _accelerator.Allocate1D(1); + _firstHitBuffer = _accelerator.Allocate1D(1); + } + + bool allocateBuffers = _exponentsBuffer == null; + if (!allocateBuffers && _capacity < requiredCapacity) + { + _exponentsBuffer?.Dispose(); + _hitBuffer?.Dispose(); + if (_hostBuffer != null) + { + ThreadStaticPools.UlongPool.Return(_hostBuffer, clearArray: false); + _hostBuffer = null!; + } + + _exponentsBuffer = null!; + _hitBuffer = null!; + allocateBuffers = true; + } + + if (allocateBuffers) + { + int desiredCapacity = requiredCapacity > _capacity ? requiredCapacity : _capacity; + _capacity = desiredCapacity; + _exponentsBuffer = _accelerator.Allocate1D(_capacity); + _hitBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ThreadStaticPools.UlongPool.Rent(_capacity); + } } public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits) @@ -731,7 +557,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } int length = primes.Length; - // EvenPerfectBitScanner always supplies at least one exponent per divisor check, so the previous zero-length guard stays commented out. + // EvenPerfectBitScanner always supplies at least one exponent per divisor check, so the guard stays commented out. // if (length == 0) // { // return; @@ -745,43 +571,79 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // return; // } + ulong cycle = divisorCycle; + ulong firstPrime = primes[0]; + + if (cycle == 0UL) + { + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(divisor, firstPrime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL) + { + cycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + cycle = computedCycle; + } + + if (cycle == 0UL) + { + throw new InvalidOperationException("GPU divisor cycle length must be non-zero."); + } + } + + ulong firstRemainder = cycle == 0UL ? 0UL : firstPrime % cycle; + Monitor.Enter(_lease.ExecutionLock); - if (_kernel == null || _hitConversionKernel == null) - { - EnsureExecutionResourcesLocked(); - } + EnsureExecutionResourcesLocked(length); - var kernel = _kernel!; - var hitConversionKernel = _hitConversionKernel!; - - ArrayView1D exponentsView = _exponentsBuffer.View; - ArrayView1D resultsView = _resultsBuffer.View; + ArrayView1D divisorView = _divisorBuffer.View; + ArrayView1D offsetView = _offsetBuffer.View; + ArrayView1D countView = _countBuffer.View; + ArrayView1D cycleView = _cycleBuffer.View; + ArrayView1D remainderView = _firstRemainderBuffer.View; + ArrayView1D firstHitView = _firstHitBuffer.View; + ArrayView1D exponentView = _exponentsBuffer.View; ArrayView1D hitView = _hitBuffer.View; - Span hostSpan = _hostBuffer.AsSpan(0, _capacity); - int offset = 0; - while (offset < length) + GpuDivisorPartialData partialData = new GpuDivisorPartialData(divisor); + divisorView.SubView(0, 1).CopyFromCPU(ref partialData, 1); + + int offsetValue = 0; + offsetView.SubView(0, 1).CopyFromCPU(ref offsetValue, 1); + + int countValue = length; + countView.SubView(0, 1).CopyFromCPU(ref countValue, 1); + + ulong cycleValue = cycle; + cycleView.SubView(0, 1).CopyFromCPU(ref cycleValue, 1); + + ulong remainderValue = firstRemainder; + remainderView.SubView(0, 1).CopyFromCPU(ref remainderValue, 1); + + Span hostSpan = _hostBuffer.AsSpan(0, length); + primes.CopyTo(hostSpan); + + ArrayView1D exponentSlice = exponentView.SubView(0, length); + exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length); + + ArrayView1D hitSliceView = hitView.SubView(0, length); + + int sentinel = int.MaxValue; + firstHitView.SubView(0, 1).CopyFromCPU(ref sentinel, 1); + + var kernel = _kernel; + if (kernel is null) { - int batchSize = Math.Min(_capacity, length - offset); - ReadOnlySpan batch = primes.Slice(offset, batchSize); - batch.CopyTo(hostSpan); - - ArrayView1D exponentSlice = exponentsView.SubView(0, batchSize); - exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), batchSize); - - ArrayView1D resultSlice = resultsView.SubView(0, batchSize); - kernel(batchSize, divisorData, exponentSlice, resultSlice); - - ArrayView1D hitSliceView = hitView.SubView(0, batchSize); - hitConversionKernel(batchSize, resultSlice, hitSliceView); - - Span hitSlice = hits.Slice(offset, batchSize); - hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), batchSize); - - offset += batchSize; + Monitor.Exit(_lease.ExecutionLock); + throw new InvalidOperationException("GPU divisor kernel was not initialized."); } + kernel(1, divisorView, offsetView, countView, exponentSlice, cycleView, remainderView, hitSliceView, firstHitView); + + Span hitSlice = hits.Slice(0, length); + hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), length); + Monitor.Exit(_lease.ExecutionLock); } @@ -797,6 +659,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } + public ulong DivisorLimit { get @@ -836,50 +699,49 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return resources; } - return new BatchResources(this, accelerator, capacity); + return new BatchResources(accelerator, capacity); } private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) => _resourcePools.GetOrAdd(accelerator, static _ => []).Add(resources); private sealed class BatchResources : IDisposable { - internal BatchResources(MersenneNumberDivisorByDivisorGpuTester owner, Accelerator accelerator, int capacity) + internal BatchResources(Accelerator accelerator, int capacity) { - int cycleCapacity = DivisorCycleCache.Shared.PreferredBatchSize; - if (cycleCapacity <= 0) - { - cycleCapacity = 1; - } + int actualCapacity = Math.Max(1, capacity); - int actualCapacity = Math.Max(capacity, cycleCapacity); - - DivisorsBuffer = accelerator.Allocate1D(actualCapacity); + DivisorDataBuffer = accelerator.Allocate1D(actualCapacity); + OffsetBuffer = accelerator.Allocate1D(actualCapacity); + CountBuffer = accelerator.Allocate1D(actualCapacity); + CycleBuffer = accelerator.Allocate1D(actualCapacity); + FirstRemainderBuffer = accelerator.Allocate1D(actualCapacity); ExponentBuffer = accelerator.Allocate1D(actualCapacity); HitsBuffer = accelerator.Allocate1D(actualCapacity); HitIndexBuffer = accelerator.Allocate1D(1); - CandidateBuffer = accelerator.Allocate1D(actualCapacity); - NextStartBuffer = accelerator.Allocate1D(1); - DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); - RemainderDeltaBuffer = accelerator.Allocate1D(actualCapacity); - RemainderBuffer = accelerator.Allocate1D(actualCapacity * RemainderSlotCount); - FilterCountBuffer = accelerator.Allocate1D(1); + ArrayPool ulongPool = ThreadStaticPools.UlongPool; + ArrayPool intPool = ThreadStaticPools.IntPool; + ArrayPool partialPool = ThreadStaticPools.GpuDivisorPartialDataPool; + Divisors = ulongPool.Rent(actualCapacity); Exponents = ulongPool.Rent(actualCapacity); - DivisorData = ThreadStaticPools.MontgomeryDivisorDataPool.Rent(actualCapacity); - CycleCandidates = ulongPool.Rent(cycleCapacity); - CycleLengths = ulongPool.Rent(cycleCapacity); - RemainderDeltaKernel = owner.GetRemainderDeltaKernel(accelerator); - RemainderScanKernel = owner.GetRemainderScanKernel(accelerator); - CandidateMaskKernel = owner.GetCandidateMaskKernel(accelerator); - CandidateCompactionKernel = owner.GetCandidateCompactionKernel(accelerator); - CandidateGenerationKernel = owner.GetCandidateGenerationKernel(accelerator); - HitDetectionKernel = owner.GetHitDetectionKernel(accelerator); + DivisorData = partialPool.Rent(actualCapacity); + Offsets = intPool.Rent(actualCapacity); + Counts = intPool.Rent(actualCapacity); + Cycles = ulongPool.Rent(actualCapacity); + FirstRemainders = ulongPool.Rent(actualCapacity); Capacity = actualCapacity; - CycleCapacity = cycleCapacity; } - internal readonly MemoryBuffer1D DivisorsBuffer; + internal readonly MemoryBuffer1D DivisorDataBuffer; + + internal readonly MemoryBuffer1D OffsetBuffer; + + internal readonly MemoryBuffer1D CountBuffer; + + internal readonly MemoryBuffer1D CycleBuffer; + + internal readonly MemoryBuffer1D FirstRemainderBuffer; internal readonly MemoryBuffer1D ExponentBuffer; @@ -887,61 +749,47 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MemoryBuffer1D HitIndexBuffer; - internal readonly MemoryBuffer1D CandidateBuffer; - - internal readonly MemoryBuffer1D NextStartBuffer; - - internal readonly MemoryBuffer1D DivisorDeltaBuffer; - - internal readonly MemoryBuffer1D RemainderDeltaBuffer; - - internal readonly MemoryBuffer1D RemainderBuffer; - internal readonly MemoryBuffer1D FilterCountBuffer; - internal readonly ulong[] Divisors; internal readonly ulong[] Exponents; - internal readonly MontgomeryDivisorData[] DivisorData; + internal readonly GpuDivisorPartialData[] DivisorData; - internal readonly ulong[] CycleCandidates; + internal readonly int[] Offsets; - internal readonly ulong[] CycleLengths; + internal readonly int[] Counts; - internal readonly Action, ArrayView, byte> RemainderDeltaKernel; + internal readonly ulong[] Cycles; - internal readonly Action, ArrayView, byte, byte> RemainderScanKernel; - internal readonly Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> CandidateMaskKernel; - internal readonly Action, ArrayView, ArrayView, ArrayView> CandidateCompactionKernel; - internal readonly Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; - internal readonly Action, ArrayView> HitDetectionKernel; + internal readonly ulong[] FirstRemainders; internal readonly int Capacity; - internal readonly int CycleCapacity; - public void Dispose() { - DivisorsBuffer.Dispose(); + DivisorDataBuffer.Dispose(); + OffsetBuffer.Dispose(); + CountBuffer.Dispose(); + CycleBuffer.Dispose(); + FirstRemainderBuffer.Dispose(); ExponentBuffer.Dispose(); HitsBuffer.Dispose(); HitIndexBuffer.Dispose(); - CandidateBuffer.Dispose(); - NextStartBuffer.Dispose(); - DivisorDeltaBuffer.Dispose(); - RemainderDeltaBuffer.Dispose(); - RemainderBuffer.Dispose(); - FilterCountBuffer.Dispose(); ArrayPool ulongPool = ThreadStaticPools.UlongPool; - ArrayPool divisorPool = ThreadStaticPools.MontgomeryDivisorDataPool; + ArrayPool intPool = ThreadStaticPools.IntPool; + ArrayPool partialPool = ThreadStaticPools.GpuDivisorPartialDataPool; ulongPool.Return(Divisors, clearArray: false); ulongPool.Return(Exponents, clearArray: false); - divisorPool.Return(DivisorData, clearArray: false); - ulongPool.Return(CycleCandidates, clearArray: false); - ulongPool.Return(CycleLengths, clearArray: false); + partialPool.Return(DivisorData, clearArray: false); + intPool.Return(Offsets, clearArray: false); + intPool.Return(Counts, clearArray: false); + ulongPool.Return(Cycles, clearArray: false); + ulongPool.Return(FirstRemainders, clearArray: false); } } + + private sealed class AcceleratorReferenceComparer : IEqualityComparer { public static AcceleratorReferenceComparer Instance { get; } = new(); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index c3fa2b1..772fab8 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -81,7 +81,7 @@ public class MersenneDivisorCycles return heuristicCycle == exponent; } - if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL) + if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL) { return cycleLength == exponent; } @@ -91,10 +91,12 @@ public class MersenneDivisorCycles public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) { - if (!IsValidMersenneDivisorCandidate(divisor, exponent)) - { - return false; - } + // EvenPerfectBitScanner only calls this helper with divisors generated from 2 * k * p + 1 and primes greater than one, + // so the validation guard would never trigger in production. Keep it commented out to leave the hot path branch-free. + // if (!IsValidMersenneDivisorCandidate(divisor, exponent)) + // { + // return false; + // } if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) { @@ -245,7 +247,7 @@ public class MersenneDivisorCycles return order; } - public static bool TryCalculateCycleLengthForExponent( + public static bool TryCalculateCycleLengthForExponentCpu( ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, @@ -255,17 +257,19 @@ public class MersenneDivisorCycles cycleLength = 0UL; primeOrderFailed = false; - if ((divisor & (divisor - 1UL)) == 0UL) - { - cycleLength = 1UL; - return true; - } + // EvenPerfectBitScanner never feeds powers of two into this helper, so leave the branch commented out to avoid redundant work. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // cycleLength = 1UL; + // return true; + // } - if (divisor <= 3UL || exponent <= 1UL) - { - cycleLength = CalculateCycleLength(divisor, divisorData); - return true; - } + // The scanner only calls this path with prime exponents greater than one and odd divisors above three, keeping the small-parameter guard unreachable. + // if (divisor <= 3UL || exponent <= 1UL) + // { + // cycleLength = CalculateCycleLength(divisor, divisorData); + // return true; + // } ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, @@ -482,7 +486,7 @@ public class MersenneDivisorCycles ulong[]? heapCandidateArray = null; bool[]? heapEvaluationArray = null; - ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData); + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); ulong order = initialOrder; for (int i = 0; i < count; i++) @@ -513,7 +517,7 @@ public class MersenneDivisorCycles ProcessReduceOrderPrime(heapCandidateArray.AsSpan(0, multiplicity), heapEvaluationArray.AsSpan(0, multiplicity), ref order, prime, multiplicity, ref stepper); } - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); ThreadStaticPools.UlongPool.Return(primes); if (heapCandidateArray is not null) @@ -525,7 +529,7 @@ public class MersenneDivisorCycles return order; } - private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepper stepper) + private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper) { ulong working = order; int actual = 0; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 5372881..ac93691 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -163,12 +163,12 @@ internal static partial class PrimeOrderCalculator Span candidates = buffer[..actual]; candidates.Sort(); - ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData); + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); int candidateCount = candidates.Length; if (stepper.InitializeCpuIsUnity(candidates[0])) { - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); return false; } @@ -176,12 +176,12 @@ internal static partial class PrimeOrderCalculator { if (stepper.ComputeNextIsUnity(candidates[i])) { - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); return false; } } - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); return true; } @@ -219,7 +219,7 @@ internal static partial class PrimeOrderCalculator buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value)); - ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData); + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); const int StackExponentCapacity = 16; const int ExponentHardLimit = 256; @@ -256,7 +256,7 @@ internal static partial class PrimeOrderCalculator ProcessExponentLoweringPrime(heapCandidateArray.AsSpan(0, exponent), heapEvaluationArray.AsSpan(0, exponent), ref order, primeFactor, exponent, ref stepper); } - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); if (heapCandidateArray is not null) { @@ -268,7 +268,7 @@ internal static partial class PrimeOrderCalculator return order; } - private static void ProcessExponentLoweringPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepper stepper) + private static void ProcessExponentLoweringPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepperCpu stepper) { ulong working = order; int actual = 0; @@ -388,7 +388,7 @@ internal static partial class PrimeOrderCalculator Span stackBuffer = stackalloc ulong[StackExponentCapacity]; ulong[]? heapCandidateArray = null; - ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData); + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); bool violates = false; for (int i = 0; i < length; i++) @@ -425,7 +425,7 @@ internal static partial class PrimeOrderCalculator } } - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); if (heapCandidateArray is not null) { @@ -435,7 +435,7 @@ internal static partial class PrimeOrderCalculator return !violates; } - private static bool ValidateOrderForFactor(Span buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepper stepper) + private static bool ValidateOrderForFactor(Span buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepperCpu stepper) { ulong working = order; int actual = 0; @@ -547,7 +547,7 @@ internal static partial class PrimeOrderCalculator int candidateCount = candidates.Count; bool allowGpuBatch = true; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepper powStepper = ThreadStaticPools.RentExponentStepper(divisorData); + ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); bool powStepperInitialized = false; // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); @@ -626,7 +626,7 @@ internal static partial class PrimeOrderCalculator candidates.Clear(); ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepper(powStepper); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); result = candidate; return true; } @@ -666,7 +666,7 @@ internal static partial class PrimeOrderCalculator candidates.Clear(); ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepper(powStepper); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); result = candidate; return true; } @@ -681,7 +681,7 @@ internal static partial class PrimeOrderCalculator candidates.Clear(); ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepper(powStepper); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); // DebugLog("No candidate confirmed"); return false; } @@ -892,7 +892,7 @@ internal static partial class PrimeOrderCalculator Span stackBuffer = stackalloc ulong[StackExponentCapacity]; ulong[]? heapCandidateArray = null; - ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData); + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); bool violates = false; for (int i = 0; i < length; i++) @@ -929,7 +929,7 @@ internal static partial class PrimeOrderCalculator } } - ThreadStaticPools.ReturnExponentStepper(stepper); + ThreadStaticPools.ReturnExponentStepperCpu(stepper); if (heapCandidateArray is not null) { @@ -944,7 +944,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool CheckCandidateViolation(Span buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepper stepper) + private static bool CheckCandidateViolation(Span buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepperCpu stepper) { ulong working = candidate; int actual = 0; diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 2da421e..b019cf8 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -2,6 +2,7 @@ using System.Buffers; using System.Collections.Generic; using System.Runtime.CompilerServices; using PerfectNumbers.Core.Cpu; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core { @@ -102,6 +103,18 @@ namespace PerfectNumbers.Core } } + [ThreadStatic] + private static ArrayPool? _gpuDivisorPartialDataPool; + + internal static ArrayPool GpuDivisorPartialDataPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _gpuDivisorPartialDataPool ??= ArrayPool.Create(); + } + } + [ThreadStatic] private static ArrayPool? _intPool; @@ -329,31 +342,31 @@ namespace PerfectNumbers.Core } [ThreadStatic] - private static ExponentRemainderStepper _exponentRemainderStepper; + private static ExponentRemainderStepperCpu _exponentRemainderStepperCpu; [ThreadStatic] - private static bool _hasExponentRemainderStepper; + private static bool _hasExponentRemainderStepperCpu; [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ExponentRemainderStepper RentExponentStepper(in MontgomeryDivisorData divisorData) + internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData) { - if (_hasExponentRemainderStepper && _exponentRemainderStepper.MatchesDivisor(divisorData)) + if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData)) { - ExponentRemainderStepper stepper = _exponentRemainderStepper; - _hasExponentRemainderStepper = false; + ExponentRemainderStepperCpu stepper = _exponentRemainderStepperCpu; + _hasExponentRemainderStepperCpu = false; stepper.Reset(); return stepper; } - return new ExponentRemainderStepper(divisorData); + return new ExponentRemainderStepperCpu(divisorData); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static void ReturnExponentStepper(ExponentRemainderStepper stepper) + internal static void ReturnExponentStepperCpu(ExponentRemainderStepperCpu stepper) { stepper.Reset(); - _exponentRemainderStepper = stepper; - _hasExponentRemainderStepper = true; + _exponentRemainderStepperCpu = stepper; + _hasExponentRemainderStepperCpu = true; } } } diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index d831097..6b736ee 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -1,3 +1,4 @@ +using System; using System.Collections.Concurrent; using System.Runtime.CompilerServices; using ILGPU; @@ -18,6 +19,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] + [Obsolete("Use ULongExtensions.MulModGpu for GPU-compatible host code or GpuUInt128.MulMod inside kernels.")] public static ulong MulMod64Gpu(this ulong a, ulong b, ulong modulus) { // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, @@ -27,6 +29,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] + [Obsolete("Use ULongExtensions.MulModGpu for GPU-compatible host code or GpuUInt128.MulMod inside kernels.")] public static ulong MulMod64GpuDeferred(this ulong a, ulong b, ulong modulus) { // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the @@ -35,6 +38,13 @@ public static partial class ULongExtensions return state.MulModWithNativeModulo(b, modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulModGpu(this ulong left, ulong right, ulong modulus) + { + GpuUInt128 state = new(left); + return state.MulMod(right, modulus); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) { @@ -71,7 +81,7 @@ public static partial class ULongExtensions while (index >= 0) { ulong currentBit = (exponent >> index) & 1UL; - ulong squared = result.MulMod64Gpu(result, modulus); + ulong squared = result.MulModGpu(result, modulus); result = currentBit == 0UL ? squared : result; index = currentBit == 0UL ? index - 1 : index; if (currentBit == 0UL) @@ -91,14 +101,14 @@ public static partial class ULongExtensions int windowLength = index - windowStart + 1; for (int square = 0; square < windowLength; square++) { - result = result.MulMod64Gpu(result, modulus); + result = result.MulModGpu(result, modulus); } ulong mask = (1UL << windowLength) - 1UL; ulong windowValue = (exponent >> windowStart) & mask; int tableIndex = (int)((windowValue - 1UL) >> 1); ulong multiplier = oddPowers[tableIndex]; - result = result.MulMod64Gpu(multiplier, modulus); + result = result.MulModGpu(multiplier, modulus); index = windowStart - 1; } @@ -108,6 +118,79 @@ public static partial class ULongExtensions return result; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong Pow2ModWindowedGpuKernel(ulong exponent, ulong modulus) + { + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + ReadOnlyGpuUInt128 modulusWide = new(modulus); + GpuUInt128 result = new(1UL); + int index = bitLength - 1; + + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + if (currentBit == 0UL) + { + result.MulMod(result.AsReadOnly(), modulusWide); + index--; + continue; + } + + int windowStartCandidate = index - windowSize + 1; + int negativeMask = windowStartCandidate >> 31; + windowStartCandidate &= ~negativeMask; + int windowStart = GetNextSetBitIndexGpu(exponent, windowStartCandidate); + int windowLength = index - windowStart + 1; + + for (int square = 0; square < windowLength; square++) + { + result.MulMod(result.AsReadOnly(), modulusWide); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + ulong multiplier = ComputeWindowedOddPowerGpuKernel(windowValue, modulusWide); + result.MulMod(multiplier, modulusWide); + + index = windowStart - 1; + } + + return result.Low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeWindowedOddPowerGpuKernel(ulong windowValue, in ReadOnlyGpuUInt128 modulus) + { + GpuUInt128 result = new(2UL); + if (windowValue == 1UL) + { + // Single-bit windows show up whenever the bit scan finds an isolated exponent bit, + // which happens routinely on the GPU path, so keep the fast return enabled. + return 2UL; + } + + ulong remaining = (windowValue - 1UL) >> 1; + GpuUInt128 squareBase = new(2UL); + squareBase.MulMod(squareBase.AsReadOnly(), modulus); + + while (remaining != 0UL) + { + if ((remaining & 1UL) != 0UL) + { + result.MulMod(squareBase.AsReadOnly(), modulus); + } + + remaining >>= 1; + if (remaining != 0UL) + { + squareBase.MulMod(squareBase.AsReadOnly(), modulus); + } + } + + return result.Low; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorData divisor, ulong exponent) { @@ -228,11 +311,11 @@ public static partial class ULongExtensions return; } - ulong square = baseValue.MulMod64Gpu(baseValue, modulus); + ulong square = baseValue.MulModGpu(baseValue, modulus); for (int i = 1; i < oddPowers.Length; i++) { ulong previous = oddPowers[i - 1]; - oddPowers[i] = previous.MulMod64Gpu(square, modulus); + oddPowers[i] = previous.MulModGpu(square, modulus); } }