From 6217ea7888358b26c1525387a91c4894454db343 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 18 Oct 2025 17:42:33 +0200 Subject: [PATCH] Document GPU kernel audit progress --- .../Gpu/GpuKernelCompilationTests.cs | 21 +- PerfectNumbers.Core/DivisorCycleCache.cs | 75 +- PerfectNumbers.Core/Gpu/GpuKernelLease.cs | 528 +------ .../Gpu/KernelAuditProgress.md | 16 + .../Gpu/Kernels/DivisorByDivisorKernels.cs | 153 ++ .../Gpu/Kernels/DivisorCycleKernels.cs | 170 +++ .../Gpu/Kernels/DivisorKernels.cs | 51 + .../Gpu/Kernels/IncrementalKernels.cs | 219 +++ .../Gpu/Kernels/KernelMathHelpers.cs | 69 + .../Gpu/Kernels/LucasLehmerKernels.cs | 452 ++++++ .../Gpu/Kernels/NttBarrettKernels.cs | 226 +++ .../Gpu/Kernels/NttButterflyKernels.cs | 30 + .../Gpu/Kernels/NttMontgomeryKernels.cs | 109 ++ .../Gpu/Kernels/NttPointwiseKernels.cs | 20 + .../Gpu/Kernels/NttTransformKernels.cs | 69 + .../Gpu/Kernels/NumericProbeKernels.cs | 51 + .../Gpu/Kernels/OrderKernels.cs | 51 + .../Gpu/Kernels/Pow2ModKernels.cs | 182 +++ .../Gpu/Kernels/Pow2MontgomeryKernels.cs | 17 + .../Gpu/Kernels/PrimeOrderKernels.cs | 1358 +++++++++++++++++ .../Gpu/Kernels/PrimeTesterKernels.cs | 99 ++ .../Gpu/Kernels/SmallPrimeFactorKernels.cs | 66 + ...MersenneNumberDivisorByDivisorGpuTester.cs | 155 +- .../Gpu/MersenneNumberDivisorGpuTester.cs | 45 +- .../Gpu/MersenneNumberLucasLehmerGpuTester.cs | 468 +----- PerfectNumbers.Core/Gpu/NttGpuMath.cs | 452 +----- .../Gpu/PrimeOrderGpuHeuristics.cs | 1347 +--------------- PerfectNumbers.Core/MersenneDivisorCycles.cs | 91 +- PerfectNumbers.Core/PrimeTester.cs | 94 +- PerfectNumbers.Core/ULongExtensions.cs | 7 +- tools/GpuNumericTypeProbe.csx | 51 +- 31 files changed, 3496 insertions(+), 3246 deletions(-) create mode 100644 PerfectNumbers.Core/Gpu/KernelAuditProgress.md create mode 100644 PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/DivisorCycleKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/DivisorKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/KernelMathHelpers.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/LucasLehmerKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/NttBarrettKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/NttButterflyKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/NttMontgomeryKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/NttPointwiseKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/NttTransformKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/NumericProbeKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/OrderKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs create mode 100644 PerfectNumbers.Core/Gpu/Kernels/SmallPrimeFactorKernels.cs diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index 7729e10..1ef8592 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -131,22 +131,13 @@ public class GpuKernelCompilationTests ?? throw new InvalidOperationException("PrimeOrderGpuHeuristics.GetPow2ModWideKernel not found."); private static readonly Action, ArrayView1D> MersenneDivisorCyclesKernel = - (Action, ArrayView1D>)Delegate.CreateDelegate( - typeof(Action, ArrayView1D>), - typeof(MersenneDivisorCycles).GetMethod("GpuDivisorCycleKernel", BindingFlags.NonPublic | BindingFlags.Static) - ?? throw new InvalidOperationException("MersenneDivisorCycles.GpuDivisorCycleKernel not found.")); + DivisorCycleKernels.GpuDivisorCycleKernel; private static readonly Action, ArrayView, ArrayView> PrimeTesterSmallPrimeKernel = - (Action, ArrayView, ArrayView>)Delegate.CreateDelegate( - typeof(Action, ArrayView, ArrayView>), - typeof(PrimeTester).GetMethod("SmallPrimeSieveKernel", BindingFlags.NonPublic | BindingFlags.Static) - ?? throw new InvalidOperationException("PrimeTester.SmallPrimeSieveKernel not found.")); + PrimeTesterKernels.SmallPrimeSieveKernel; private static readonly Action, ArrayView> PrimeTesterSharesFactorKernel = - (Action, ArrayView>)Delegate.CreateDelegate( - typeof(Action, ArrayView>), - typeof(PrimeTester).GetMethod("SharesFactorKernel", BindingFlags.NonPublic | BindingFlags.Static) - ?? throw new InvalidOperationException("PrimeTester.SharesFactorKernel not found.")); + PrimeTesterKernels.SharesFactorKernel; public static IEnumerable KernelLoaders() { @@ -180,9 +171,9 @@ public class GpuKernelCompilationTests yield return Loader("PrimeOrderGpuHeuristics.CalculateOrderKernel", CompilePrimeOrderKernel); yield return Loader("PrimeOrderGpuHeuristics.Pow2ModKernel", CompilePrimeOrderPow2ModKernel); yield return Loader("PrimeOrderGpuHeuristics.Pow2ModKernelWide", CompilePrimeOrderPow2ModWideKernel); - yield return Loader("MersenneDivisorCycles.GpuDivisorCycleKernel", CompileMersenneDivisorCyclesKernel); - yield return Loader("PrimeTester.SmallPrimeSieveKernel", CompilePrimeTesterSmallPrimeKernel); - yield return Loader("PrimeTester.SharesFactorKernel", CompilePrimeTesterSharesFactorKernel); + yield return Loader("DivisorCycleKernels.GpuDivisorCycleKernel", CompileMersenneDivisorCyclesKernel); + yield return Loader("PrimeTesterKernels.SmallPrimeSieveKernel", CompilePrimeTesterSmallPrimeKernel); + yield return Loader("PrimeTesterKernels.SharesFactorKernel", CompilePrimeTesterSharesFactorKernel); yield return Loader("GpuKernelPool.KernelLeaseKernels", CompileGpuKernelPoolKernels); } diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index ba320f7..801d51e 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -269,80 +269,7 @@ public sealed class DivisorCycleCache private static Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> LoadKernel(Accelerator accelerator) { - return accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(GpuAdvanceDivisorCyclesKernel); - } - - private static void GpuAdvanceDivisorCyclesKernel( - Index1D index, - int steps, - ArrayView1D divisors, - ArrayView1D pow, - ArrayView1D order, - ArrayView1D cycles, - ArrayView1D status) - { - if (status[index] != ByteZero) - { - return; - } - - ulong divisor = divisors[index]; - - if ((divisor & (divisor - 1UL)) == 0UL) - { - cycles[index] = 1UL; - status[index] = ByteOne; - pow[index] = 1UL; - order[index] = 1UL; - return; - } - - ulong currentPow = pow[index]; - ulong currentOrder = order[index]; - - if (divisor <= 3UL) - { - while (currentPow != 1UL) - { - currentPow += currentPow; - if (currentPow >= divisor) - { - currentPow -= divisor; - } - - currentOrder++; - } - - cycles[index] = currentOrder; - status[index] = ByteOne; - pow[index] = currentPow; - order[index] = currentOrder; - return; - } - - do - { - currentPow += currentPow; - if (currentPow >= divisor) - { - currentPow -= divisor; - } - - currentOrder++; - - if (currentPow == 1UL) - { - cycles[index] = currentOrder; - status[index] = ByteOne; - pow[index] = currentPow; - order[index] = currentOrder; - return; - } - } - while (--steps != 0); - - pow[index] = currentPow; - order[index] = currentOrder; + return accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(DivisorCycleKernels.GpuAdvanceDivisorCyclesKernel); } private sealed class AcceleratorReferenceComparer : IEqualityComparer diff --git a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs index 914f026..f590795 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelLease.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelLease.cs @@ -1,5 +1,4 @@ using System.Collections.Concurrent; -using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; using static PerfectNumbers.Core.Gpu.GpuContextPool; @@ -66,7 +65,7 @@ public sealed class GpuKernelLease { return KernelContainer.InitOnce(ref _kernels!.Order, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernelScan); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernels.OrderKernelScan); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView>>(); }); @@ -83,7 +82,7 @@ public sealed class GpuKernelLease { return KernelContainer.InitOnce(ref _kernels!.Pow2Mod, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernels.Pow2ModKernelScan); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); }); @@ -100,7 +99,7 @@ public sealed class GpuKernelLease { return KernelContainer.InitOnce(ref _kernels!.Incremental, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernelScan); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernels.IncrementalKernelScan); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView1D>>(); }); @@ -117,7 +116,7 @@ public sealed class GpuKernelLease { return KernelContainer.InitOnce(ref _kernels!.IncrementalOrder, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalOrderKernelScan); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernels.IncrementalOrderKernelScan); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView1D>>(); }); @@ -134,7 +133,7 @@ public sealed class GpuKernelLease { return KernelContainer.InitOnce(ref _kernels!.Pow2ModOrder, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModOrderKernelScan); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModKernels.Pow2ModOrderKernelScan); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView1D>>(); }); @@ -151,7 +150,7 @@ public sealed class GpuKernelLease { return KernelContainer.InitOnce(ref _kernels!.SmallPrimeFactor, () => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernelScan); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); }); @@ -159,211 +158,6 @@ public sealed class GpuKernelLease } } - private static void SmallPrimeFactorKernelScan( - Index1D index, - ulong value, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D primeSlots, - ArrayView1D exponentSlots, - ArrayView1D countSlot, - ArrayView1D remainingSlot) - { - if (index != 0) - { - return; - } - - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - ulong remaining = value; - int count = 0; - int capacity = (int)primeSlots.Length; - - for (int i = 0; i < primeCount && remaining > 1UL && count < capacity; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate == 0U || primeCandidate > effectiveLimit) - { - break; - } - - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remaining) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == 0UL); - - primeSlots[count] = primeValue; - exponentSlots[count] = exponent; - count++; - } - - countSlot[0] = count; - remainingSlot[0] = remaining; - } - - // TODO: Plumb the small-cycles device buffer into all kernels that can benefit - // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. - - private static void IncrementalKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, - ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView orders, - ArrayView1D smallCycles) - { - ulong idx = (ulong)index.X; - // TODO: Replace these `%` computations with the precomputed Mod3/Mod5 tables so GPU kernels reuse cached residues - // instead of performing modulo operations that the benchmarks showed slower on wide sweeps. - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - // residue automaton - ulong step10 = (exponent.Mod10() << 1).Mod10(); - ulong r10 = q0m10 + (step10 * idx).Mod10(); - r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5 - if (shouldCheck) - { - ulong step8 = ((exponent & 7UL) << 1) & 7UL; - ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked - // bitmask-based implementation for GPU workloads. - ulong step3 = ((exponent % 3UL) << 1) % 3UL; - ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong step5Loc = ((exponent % 5UL) << 1) % 5UL; - ulong r5 = q0m5 + (step5Loc * idxMod5); if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - - if (!shouldCheck) - { - orders[index] = 0UL; - return; - } - - // TODO: Is this expected that the value of kStart is modified? - GpuUInt128 k = kStart + (GpuUInt128)idx; - GpuUInt128 q = twoP; - q.Mul64(k); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - - // Small-cycles in-kernel early rejection from device table - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - orders[index] = 0UL; - return; - } - } - - // TODO: q will be lowered by GpuUInt128.One here. Is this expected? Maybe that's why we have issues with the incremental kernel? Review other places where we use GpuUInt128 operators. - GpuUInt128 phi = q - GpuUInt128.One; - // TODO: Do we have the check here because we don't support high values or this is expected result? - if (phi.High != 0UL) - { - orders[index] = 0UL; - return; - } - - ulong phi64 = phi.Low; - // TODO: Replace these Pow2Mod calls with the ProcessEightBitWindows helper when the shared windowed - // scalar implementation lands; benchmarks showed the windowed kernel trimming per-divisor runtime by ~2.4×. - if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - ulong div = FastDiv64Gpu(phi64, exponent, divMul); - GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - orders[index] = 1UL; - } - - private static void OrderKernelScan(Index1D index, ulong exponent, ulong divMul, ArrayView qs, ArrayView orders) - { - GpuUInt128 q = qs[index]; - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - GpuUInt128 phi = q - GpuUInt128.One; - if (phi.High != 0UL) - { - orders[index] = 0UL; - return; - } - - ulong phi64 = phi.Low; - // TODO: Once the ProcessEightBitWindows helper is available, switch this order kernel to that faster - // Pow2Mod variant so cycle checks inherit the same gains observed in GpuPow2ModBenchmarks. - GpuUInt128 pow = GpuUInt128.Pow2Mod(phi64, in readOnlyQ); - if (pow != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - ulong div = FastDiv64Gpu(phi64, exponent, divMul); - GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) - { - orders[index] = 0UL; - return; - } - - orders[index] = exponent; - } - private void Dispose(bool disposing) { if (_disposed) @@ -390,316 +184,6 @@ public sealed class GpuKernelLease Pool.Enqueue(this); } - private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, - ResidueAutomatonArgs ra, ArrayView orders, - ArrayView1D smallCycles, - ArrayView1D smallPrimesLastOne, - ArrayView1D smallPrimesLastSeven, - ArrayView1D smallPrimesPow2LastOne, - ArrayView1D smallPrimesPow2LastSeven) - { - ulong idx = (ulong)index.X; - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - orders[index] = 0UL; - return; - } - - // TODO: kStart is modified after this. Is this expected? - kStart.Add(idx); - GpuUInt128 q = twoP; - q.Mul64(kStart); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - // cycle should always be initialized if we're within array limit in production code - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - orders[index] = 0UL; - return; - } - } - // TODO: Swap Pow2Minus1Mod for the eight-bit window helper once the scalar version switches; - // benchmarks show the windowed variant cuts large-divisor scans from ~51 µs to ~21 µs. - if (GpuUInt128.Pow2Minus1Mod(exponent, in readOnlyQ) != GpuUInt128.Zero) - { - orders[index] = 0UL; - return; - } - - ArrayView1D primes = lastIsSeven != 0 ? smallPrimesLastSeven : smallPrimesLastOne; - ArrayView1D primesPow2 = lastIsSeven != 0 ? smallPrimesPow2LastSeven : smallPrimesPow2LastOne; - int primesLen = (int)primes.Length; - for (int i = 0; i < primesLen; i++) - { - ulong square = primesPow2[i]; - if (new GpuUInt128(0UL, square) > q) - { - break; - } - ulong prime = primes[i]; - if (Mod128By64(q, prime) == 0UL) - { - orders[index] = 0UL; - return; - } - } - - orders[index] = exponent; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong FastDiv64Gpu(ulong value, ulong divisor, ulong mul) - { - ulong quotient = GpuUInt128.MulHigh(value, mul); - GpuUInt128 remainder = new(0UL, value); - GpuUInt128 product = new(GpuUInt128.MulHigh(quotient, divisor), quotient * divisor); - remainder.Sub(product); - - if (remainder.High != 0UL || remainder.Low >= divisor) - { - quotient++; - } - - return quotient; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong Mod128By64(GpuUInt128 value, ulong modulus) - { - ulong rem = 0UL; - ulong part = value.High; - for (int i = 0; i < 64; i++) - { - rem = (rem << 1) | (part >> 63); - if (rem >= modulus) - { - rem -= modulus; - } - part <<= 1; - } - - part = value.Low; - for (int i = 0; i < 64; i++) - { - rem = (rem << 1) | (part >> 63); - if (rem >= modulus) - { - rem -= modulus; - } - part <<= 1; - } - - return rem; - } - - // Device-friendly small cycle length for 2 mod divisor - private static ulong CalculateCycleLengthSmall(ulong divisor) - { - if ((divisor & (divisor - 1UL)) == 0UL) - return 1UL; - - ulong order = 1UL, pow = 2UL; - while (pow != 1UL) - { - pow <<= 1; - if (pow > divisor) - pow -= divisor; - - order++; - } - - return order; - } - - private static void IncrementalOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, - ResidueAutomatonArgs ra, ArrayView found, ArrayView1D smallCycles) - { - ulong idx = (ulong)index.X; - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - return; - } - - GpuUInt128 k = kStart + (GpuUInt128)idx; - GpuUInt128 q = twoP; - q.Mul64(k); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - - // Small-cycles in-kernel early rejection from device table - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - return; - } - } - GpuUInt128 phi = q - GpuUInt128.One; - if (phi.High != 0UL) - { - return; - } - - ulong phi64 = phi.Low; - // TODO: Upgrade this pow2mod order kernel to the ProcessEightBitWindows helper once available so GPU residue - // scans avoid the single-bit ladder that benchmarks found to be 2.3× slower on large exponents. - if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) - { - return; - } - - GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) - { - return; - } - - ulong div = FastDiv64Gpu(phi64, exponent, divMul); - GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; - if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) - { - return; - } - - Atomic.Or(ref found[0], 1); - } - - private static void Pow2ModOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, - ResidueAutomatonArgs ra, ArrayView found, ArrayView1D smallCycles) - { - ulong idx = (ulong)index.X; - ulong idxMod3 = idx % 3UL; - ulong idxMod5 = idx % 5UL; - ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - return; - } - - GpuUInt128 k = kStart + (GpuUInt128)idx; - GpuUInt128 q = twoP; - q.Mul64(k); - q.Add(GpuUInt128.One); - ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); - // Small-cycles in-kernel early rejection from device table - if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) - { - ulong cycle = smallCycles[(int)q.Low]; - if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) - { - return; - } - } - // TODO: Replace this Pow2Mod check with the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; - // residue order scans will then benefit from the same 2× speedup the benchmarked windowed kernel delivered. - if (GpuUInt128.Pow2Mod(exponent, in readOnlyQ) != GpuUInt128.One) - { - return; - } - - Atomic.Or(ref found[0], 1); - } - public void Dispose() { Dispose(disposing: false); diff --git a/PerfectNumbers.Core/Gpu/KernelAuditProgress.md b/PerfectNumbers.Core/Gpu/KernelAuditProgress.md new file mode 100644 index 0000000..010c1dc --- /dev/null +++ b/PerfectNumbers.Core/Gpu/KernelAuditProgress.md @@ -0,0 +1,16 @@ +# GPU Kernel Refactor Audit Progress + +## Review Summary +- [x] MersenneNumberDivisorByDivisorGpuTester: Confirmed kernel bindings, resource pooling, and filtering comments (including Montgomery data handling notes) remain intact with no functional changes. +- [x] MersenneNumberDivisorGpuTester: Verified the divisibility helper, candidate bootstrap, and divisor-cycle integration preserve their original comments and TODO guidance verbatim. +- [x] MersenneNumberIncrementalGpuTester: Checked that the incremental scan workflow, residue preparation comments, and TODO placeholders stayed untouched. +- [x] MersenneNumberLucasLehmerGpuTester: Audited the Lucas–Lehmer kernels, TODO annotations, and GPU limiter usage to ensure they were moved verbatim. +- [x] MersenneNumberOrderGpuTester: Confirmed residue automaton notes, TODO entries, and kernel selections remain unchanged. +- [x] MersenneNumberResidueGpuTester: Ensured the residue scan notes about cycle integration and buffer zeroing are still present without logic edits. +- [x] NttGpuMath: Reviewed cached kernel launchers and TODO lists to verify that only structural moves occurred. +- [x] Pow2OddPowerTable (PrimeOrderGpuHeuristics partial): Checked the odd-power ladder implementation and comments for exact retention. +- [x] PrimeOrderGpuHeuristics: Validated that partial factor kernels, capability overrides, and extensive TODO documentation remain verbatim. +- [x] ReadOnlyGpuUInt128: Confirmed helper constructors and implicit conversions were untouched by the refactor. + +## Next Step +Audit complete for the GPU namespace after the kernel extraction. No further classes remain pending in this pass. diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs new file mode 100644 index 0000000..f3c0a17 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -0,0 +1,153 @@ +using ILGPU; +using ILGPU.Runtime; + +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; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits[index] = 0; + return; + } + + ulong exponent = exponents[index]; + hits[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; + } + + 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; + 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) + { + int addend = 0; + if (localIndex >= offset) + { + addend = shared[localIndex - offset]; + } + + Group.Barrier(); + + if (localIndex >= offset) + { + int sum = shared[localIndex] + addend; + if (sum >= mod) + { + sum %= mod; + } + + shared[localIndex] = sum; + } + + Group.Barrier(); + + offset <<= 1; + } + + if (globalIndex == 0) + { + remainders[0] = baseRemainder; + return; + } + + int remainder = baseRemainder + shared[localIndex]; + if (remainder >= mod) + { + remainder %= 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; + int length = (int)mask.Length; + if (globalIndex >= length) + { + return; + } + + byte value10 = remainder10[globalIndex]; + bool accept10; + if (lastIsSevenFlag != 0) + { + accept10 = value10 == 3 || value10 == 7 || value10 == 9; + } + else + { + accept10 = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9; + } + + if (!accept10) + { + mask[globalIndex] = 0; + return; + } + + byte value8 = remainder8[globalIndex]; + if (value8 != 1 && value8 != 7) + { + mask[globalIndex] = 0; + return; + } + + if (remainder3[globalIndex] == 0 || remainder5[globalIndex] == 0) + { + mask[globalIndex] = 0; + return; + } + + mask[globalIndex] = 1; + } + + public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) + { + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + results[index] = 0UL; + return; + } + + ulong exponent = exponents[index]; + + results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true); + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorCycleKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorCycleKernels.cs new file mode 100644 index 0000000..e0d13ca --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorCycleKernels.cs @@ -0,0 +1,170 @@ +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class DivisorCycleKernels +{ + private const byte ByteZero = 0; + private const byte ByteOne = 1; + + // GPU kernel for divisor cycle calculation + public static void GpuDivisorCycleKernel( + Index1D index, + ArrayView1D divisors, + ArrayView1D outCycles) + { + int i = index.X; + outCycles[i] = CalculateCycleLengthGpu(divisors[i]); + } + + // GPU-friendly version of cycle length calculation + /// + /// GPU-friendly cycle calculator that unrolls sixteen doubling steps; it wins the 8,388,607 benchmark and stays within ~2% + /// of the octo loop at divisor 131,071. + /// + public static ulong CalculateCycleLengthGpu(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + return 1UL; + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool GpuStep(ref ulong pow, ulong divisor, ref ulong order) + { + pow += pow; + if (pow >= divisor) + pow -= divisor; + + order++; + return pow == 1UL; + } + + public static void GpuAdvanceDivisorCyclesKernel( + Index1D index, + int steps, + ArrayView1D divisors, + ArrayView1D pow, + ArrayView1D order, + ArrayView1D cycles, + ArrayView1D status) + { + if (status[index] != ByteZero) + { + return; + } + + ulong divisor = divisors[index]; + + if ((divisor & (divisor - 1UL)) == 0UL) + { + cycles[index] = 1UL; + status[index] = ByteOne; + pow[index] = 1UL; + order[index] = 1UL; + return; + } + + ulong currentPow = pow[index]; + ulong currentOrder = order[index]; + + if (divisor <= 3UL) + { + while (currentPow != 1UL) + { + currentPow += currentPow; + if (currentPow >= divisor) + { + currentPow -= divisor; + } + + currentOrder++; + } + + cycles[index] = currentOrder; + status[index] = ByteOne; + pow[index] = currentPow; + order[index] = currentOrder; + return; + } + + do + { + currentPow += currentPow; + if (currentPow >= divisor) + { + currentPow -= divisor; + } + + currentOrder++; + + if (currentPow == 1UL) + { + cycles[index] = currentOrder; + status[index] = ByteOne; + pow[index] = currentPow; + order[index] = currentOrder; + return; + } + } + while (--steps != 0); + + pow[index] = currentPow; + order[index] = currentOrder; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorKernels.cs new file mode 100644 index 0000000..d719c0b --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorKernels.cs @@ -0,0 +1,51 @@ +using ILGPU; +using ILGPU.Algorithms; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class DivisorKernels +{ + public static void Kernel(Index1D _, ulong exponent, ReadOnlyGpuUInt128 divisor, ArrayView result) + { + if (divisor.High == 0UL && divisor.Low <= 1UL) + { + result[0] = 0; + return; + } + + GpuUInt128 baseVal; + int x = divisor.High == 0UL + ? 63 - XMath.LeadingZeroCount(divisor.Low) + : 127 - XMath.LeadingZeroCount(divisor.High); + ReadOnlyGpuUInt128 readOnlyMod = divisor; + if (x == 0) + { + // TODO: Swap Pow2Mod for the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; + // GpuPow2ModBenchmarks showed the windowed kernel cutting large divisors from ~51 µs to ~21 µs. + baseVal = GpuUInt128.Pow2Mod(exponent, in readOnlyMod); + result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0; + return; + } + + ulong ux = (ulong)x; + if (x >= 64) + { + baseVal = new GpuUInt128(1UL << (x - 64), 0UL); + } + else + { + baseVal = new GpuUInt128(0UL, 1UL << x); + } + + GpuUInt128 pow = baseVal; + pow.ModPow(exponent / ux, in readOnlyMod); + // TODO: Replace this trailing Pow2Mod call with the ProcessEightBitWindows variant once the shared + // windowed helper lands so mixed-radix decompositions benefit from the same GPU speedup. + var part2 = GpuUInt128.Pow2Mod(exponent % ux, in readOnlyMod); + GpuUInt128 product = pow; + ReadOnlyGpuUInt128 part2ReadOnly = part2.AsReadOnly(); + product.MulMod(in part2ReadOnly, in readOnlyMod); + result[0] = product.High == 0UL && product.Low == 1UL ? (byte)1 : (byte)0; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs new file mode 100644 index 0000000..c4d0245 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs @@ -0,0 +1,219 @@ +using ILGPU; +using ILGPU.Algorithms; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class IncrementalKernels +{ + // TODO: Plumb the small-cycles device buffer into all kernels that can benefit + // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. + public static void IncrementalKernelScan( + Index1D index, + ulong exponent, + GpuUInt128 twoP, + GpuUInt128 kStart, + byte lastIsSeven, + ulong divMul, + ulong q0m10, + ulong q0m8, + ulong q0m3, + ulong q0m5, + ArrayView orders, + ArrayView1D smallCycles) + { + ulong idx = (ulong)index.X; + // TODO: Replace these `%` computations with the precomputed Mod3/Mod5 tables so GPU kernels reuse cached residues + // instead of performing modulo operations that the benchmarks showed slower on wide sweeps. + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + // residue automaton + ulong step10 = (exponent.Mod10() << 1).Mod10(); + ulong r10 = q0m10 + (step10 * idx).Mod10(); + r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5 + if (shouldCheck) + { + ulong step8 = ((exponent & 7UL) << 1) & 7UL; + ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked + // bitmask-based implementation for GPU workloads. + ulong step3 = ((exponent % 3UL) << 1) % 3UL; + ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong step5Loc = ((exponent % 5UL) << 1) % 5UL; + ulong r5 = q0m5 + (step5Loc * idxMod5); if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + + if (!shouldCheck) + { + orders[index] = 0UL; + return; + } + + // TODO: Is this expected that the value of kStart is modified? + GpuUInt128 k = kStart + (GpuUInt128)idx; + GpuUInt128 q = twoP; + q.Mul64(k); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + + // Small-cycles in-kernel early rejection from device table + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + orders[index] = 0UL; + return; + } + } + + // TODO: q will be lowered by GpuUInt128.One here. Is this expected? Maybe that's why we have issues with the incremental kernel? Review other places where we use GpuUInt128 operators. + GpuUInt128 phi = q - GpuUInt128.One; + // TODO: Do we have the check here because we don't support high values or this is expected result? + if (phi.High != 0UL) + { + orders[index] = 0UL; + return; + } + + ulong phi64 = phi.Low; + // TODO: Replace these Pow2Mod calls with the ProcessEightBitWindows helper when the shared windowed + // scalar implementation lands; benchmarks showed the windowed kernel trimming per-divisor runtime by ~2.4×. + if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + ulong div = KernelMathHelpers.FastDiv64Gpu(phi64, exponent, divMul); + GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + orders[index] = 1UL; + } + + public static void IncrementalOrderKernelScan( + Index1D index, + ulong exponent, + GpuUInt128 twoP, + GpuUInt128 kStart, + byte lastIsSeven, + ulong divMul, + ResidueAutomatonArgs ra, + ArrayView found, + ArrayView1D smallCycles) + { + ulong idx = (ulong)index.X; + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + return; + } + + GpuUInt128 k = kStart + (GpuUInt128)idx; + GpuUInt128 q = twoP; + q.Mul64(k); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + + // Small-cycles in-kernel early rejection from device table + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + return; + } + } + GpuUInt128 phi = q - GpuUInt128.One; + if (phi.High != 0UL) + { + return; + } + + ulong phi64 = phi.Low; + // TODO: Upgrade this pow2mod order kernel to the ProcessEightBitWindows helper once available so GPU residue + // scans avoid the single-bit ladder that benchmarks found to be 2.3× slower on large exponents. + if (GpuUInt128.Pow2Mod(phi64, in readOnlyQ) != GpuUInt128.One) + { + return; + } + + GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) + { + return; + } + + ulong div = KernelMathHelpers.FastDiv64Gpu(phi64, exponent, divMul); + GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) + { + return; + } + + Atomic.Or(ref found[0], 1); + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/KernelMathHelpers.cs b/PerfectNumbers.Core/Gpu/Kernels/KernelMathHelpers.cs new file mode 100644 index 0000000..1cf897f --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/KernelMathHelpers.cs @@ -0,0 +1,69 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core.Gpu; + +internal static class KernelMathHelpers +{ + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong FastDiv64Gpu(ulong value, ulong divisor, ulong mul) + { + ulong quotient = GpuUInt128.MulHigh(value, mul); + GpuUInt128 remainder = new(0UL, value); + GpuUInt128 product = new(GpuUInt128.MulHigh(quotient, divisor), quotient * divisor); + remainder.Sub(product); + + if (remainder.High != 0UL || remainder.Low >= divisor) + { + quotient++; + } + + return quotient; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Mod128By64(GpuUInt128 value, ulong modulus) + { + ulong rem = 0UL; + ulong part = value.High; + for (int i = 0; i < 64; i++) + { + rem = (rem << 1) | (part >> 63); + if (rem >= modulus) + { + rem -= modulus; + } + part <<= 1; + } + + part = value.Low; + for (int i = 0; i < 64; i++) + { + rem = (rem << 1) | (part >> 63); + if (rem >= modulus) + { + rem -= modulus; + } + part <<= 1; + } + + return rem; + } + + public static ulong CalculateCycleLengthSmall(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + return 1UL; + + ulong order = 1UL, pow = 2UL; + while (pow != 1UL) + { + pow <<= 1; + if (pow > divisor) + pow -= divisor; + + order++; + } + + return order; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/LucasLehmerKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/LucasLehmerKernels.cs new file mode 100644 index 0000000..82cbc1c --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/LucasLehmerKernels.cs @@ -0,0 +1,452 @@ +using System; +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Runtime; +using UInt128 = System.UInt128; + +namespace PerfectNumbers.Core.Gpu; + +internal static class LucasLehmerKernels +{ + public static void AddSmallKernel(Index1D index, ulong add, ArrayView value) + { + ulong carry = add; + for (int i = 0; i < value.Length && carry != 0UL; i++) + { + var limb = value[i]; + ulong low = limb.Low + carry; + ulong carryOut = low < carry ? 1UL : 0UL; + ulong high = limb.High + carryOut; + carry = high < limb.High ? 1UL : 0UL; + value[i] = new GpuUInt128(high, low); + } + } + + public static void SubtractSmallKernel(Index1D index, ulong subtract, ArrayView value) + { + ulong borrow = subtract; + for (int i = 0; i < value.Length && borrow != 0UL; i++) + { + var limb = value[i]; + ulong low = limb.Low; + ulong high = limb.High; + if (low >= borrow) + { + low -= borrow; + borrow = 0UL; + } + else + { + low = unchecked(low - borrow); + borrow = 1UL; + if (high != 0UL) + { + high--; + borrow = 0UL; + } + else + { + high = ulong.MaxValue; + } + } + + value[i] = new GpuUInt128(high, low); + } + } + + public static void ReduceModMersenneKernel(Index1D index, ulong exponent, ArrayView value) + { + int limbCount = (int)((exponent + 127UL) / 128UL); + for (int i = limbCount; i < value.Length; i++) + { + var limb = value[i]; + if (limb.IsZero) + { + continue; + } + + int target = i - limbCount; + var sum = value[target]; + var original = sum; + sum.Add(limb); + value[target] = sum; + value[i] = new GpuUInt128(0UL, 0UL); + if (sum.CompareTo(original) < 0) + { + int j = target + 1; + while (j < limbCount) + { + var next = value[j]; + var nextOriginal = next; + next.Add(1UL); + value[j] = next; + if (next.CompareTo(nextOriginal) >= 0) + { + break; + } + + j++; + } + } + } + + int topBits = (int)exponent.Mod128(); + if (topBits == 0) + { + topBits = 128; + } + + int topIndex = limbCount - 1; + ulong topLowMask = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL; + ulong topHighMask = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL; + var top = value[topIndex]; + var carryBits = top >> topBits; + top = new(top.High & topHighMask, top.Low & topLowMask); + value[topIndex] = top; + + if (!carryBits.IsZero) + { + int j = 0; + while (j < limbCount && !carryBits.IsZero) + { + var cur = value[j]; + var curOriginal = cur; + cur.Add(carryBits); + value[j] = cur; + carryBits = cur.CompareTo(curOriginal) < 0 ? new GpuUInt128(0UL, 1UL) : new GpuUInt128(0UL, 0UL); + j++; + } + } + + bool geq = true; + for (int i = topIndex; i >= 0; i--) + { + var limb = value[i]; + ulong modHigh = i == topIndex ? topHighMask : ulong.MaxValue; + ulong modLow = i == topIndex ? topLowMask : ulong.MaxValue; + if (limb.High > modHigh || (limb.High == modHigh && limb.Low > modLow)) + { + geq = true; + break; + } + + if (limb.High < modHigh || (limb.High == modHigh && limb.Low < modLow)) + { + geq = false; + break; + } + } + + if (geq) + { + ulong borrow = 0UL; + for (int i = 0; i < limbCount - 1; i++) + { + var limb = value[i]; + ulong newLow = limb.Low - ulong.MaxValue - borrow; + borrow = limb.Low < ulong.MaxValue + borrow ? 1UL : 0UL; + ulong newHigh = limb.High - ulong.MaxValue - borrow; + borrow = limb.High < ulong.MaxValue + borrow ? 1UL : 0UL; + value[i] = new GpuUInt128(newHigh, newLow); + } + + var topLimb = value[topIndex]; + ulong newTopLow = topLimb.Low - topLowMask - borrow; + borrow = topLimb.Low < topLowMask + borrow ? 1UL : 0UL; + ulong newTopHigh = topLimb.High - topHighMask - borrow; + value[topIndex] = new GpuUInt128(newTopHigh, newTopLow); + } + + for (int i = limbCount; i < value.Length; i++) + { + value[i] = new GpuUInt128(0UL, 0UL); + } + } + + public static void IsZeroKernel(Index1D index, ArrayView value, ArrayView result) + { + byte isZero = 1; + for (int i = 0; i < value.Length; i++) + { + var limb = value[i]; + if (limb.High != 0UL || limb.Low != 0UL) + { + isZero = 0; + break; + } + } + + result[0] = isZero; + } + + private static void ReduceModMersenne(Span value, ulong exponent) + { + int limbCount = (int)((exponent + 127UL) / 128UL); + for (int i = limbCount; i < value.Length; i++) + { + var limb = value[i]; + if (limb.IsZero) + { + continue; + } + + int target = i - limbCount; + UInt128 original = value[target]; + UInt128 sum = original + (UInt128)limb; + bool carry = sum < original; + value[target] = new GpuUInt128(sum); + value[i] = new GpuUInt128(0UL, 0UL); + if (carry) + { + int j = target + 1; + while (j < limbCount) + { + UInt128 next = value[j]; + UInt128 nextSum = next + 1UL; + value[j] = new GpuUInt128(nextSum); + if (nextSum != 0UL) + { + break; + } + + j++; + } + } + } + + int topBits = (int)exponent.Mod128(); + if (topBits == 0) + { + topBits = 128; + } + + int topIndex = limbCount - 1; + ulong topLowMask = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL; + ulong topHighMask = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL; + UInt128 mask = ((UInt128)topHighMask << 64) | topLowMask; + UInt128 top = value[topIndex]; + UInt128 carryBits = top >> topBits; + top &= mask; + value[topIndex] = new GpuUInt128(top); + if (carryBits != 0UL) + { + int j = 0; + UInt128 carry = carryBits; + while (j < limbCount && carry != 0UL) + { + UInt128 cur = value[j]; + UInt128 sum = cur + carry; + bool overflow = sum < cur; + value[j] = new GpuUInt128(sum); + carry = overflow ? 1UL : 0UL; + j++; + } + } + + bool geq = true; + for (int i = topIndex; i >= 0; i--) + { + UInt128 limb = value[i]; + UInt128 modLimb = i == topIndex ? mask : UInt128.MaxValue; + if (limb > modLimb) + { + geq = true; + break; + } + + if (limb < modLimb) + { + geq = false; + break; + } + } + + if (geq) + { + ulong borrow = 0UL; + for (int i = 0; i < limbCount - 1; i++) + { + ulong low = value[i].Low; + ulong high = value[i].High; + ulong newLow = low - ulong.MaxValue - borrow; + borrow = low < ulong.MaxValue + borrow ? 1UL : 0UL; + ulong newHigh = high - ulong.MaxValue - borrow; + borrow = high < ulong.MaxValue + borrow ? 1UL : 0UL; + value[i] = new GpuUInt128(newHigh, newLow); + } + + ulong topLow = value[topIndex].Low; + ulong topHigh = value[topIndex].High; + ulong newTopLow = topLow - topLowMask - borrow; + borrow = topLow < topLowMask + borrow ? 1UL : 0UL; + ulong newTopHigh = topHigh - topHighMask - borrow; + value[topIndex] = new GpuUInt128(newTopHigh, newTopLow); + } + + for (int i = limbCount; i < value.Length; i++) + { + value[i] = new GpuUInt128(0UL, 0UL); + } + } + + public static void KernelBatch(Index1D index, ArrayView exponents, ArrayView moduli, ArrayView states) + { + int idx = index.X; + ulong exponent = exponents[idx]; + GpuUInt128 modulus = moduli[idx]; + var s = new GpuUInt128(0UL, 4UL); + ulong limit = exponent - 2UL; + for (ulong i = 0UL; i < limit; i++) + { + s = SquareModMersenne128(s, exponent); + s.SubMod(2UL, modulus); + } + + states[idx] = s; + } + + public static void Kernel(Index1D index, ulong exponent, GpuUInt128 modulus, ArrayView state) + { + // Lucas–Lehmer iteration in the field GF(2^p-1). + // For p < 128, use a Mersenne-specific squaring + reduction to avoid + // the generic 128-bit long-division-like reduction. + var s = new GpuUInt128(0UL, 4UL); + ulong i = 0UL, limit = exponent - 2UL; + if (exponent < 128UL) + { + for (; i < limit; i++) + { + s = SquareModMersenne128(s, exponent); + s.SubMod(2UL, modulus); + } + } + else + { + for (; i < limit; i++) + { + s.SquareMod(modulus); + s.SubMod(2UL, modulus); + } + } + + state[0] = s; + } + + // Squares a 128-bit value and reduces modulo M_p = 2^p - 1 using + // Mersenne folding. Valid for 0 < p < 128 and input s in [0, M_p]. + private static GpuUInt128 SquareModMersenne128(GpuUInt128 s, ulong p) + { + // Full square: (H*2^64 + L)^2 = H^2*2^128 + 2HL*2^64 + L^2 + SquareFull128(s, out var q3, out var q2, out var q1, out var q0); + + // v = low128 + high128 (fold by 128 bits once) + var vHigh = new GpuUInt128(q3, q2); + var vLow = new GpuUInt128(q1, q0); + vLow.Add(vHigh); + + // Mask to keep only p low bits and fold the carry bits (v >> p) + int topBits = (int)p.Mod128(); + if (topBits == 0) + { + topBits = 128; + } + + ulong maskLow = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL; + ulong maskHigh = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL; + + var carryBits = vLow >> topBits; + vLow = new(vLow.High & maskHigh, vLow.Low & maskLow); + + // Propagate carryBits back into the masked value (single-limb fold) + GpuUInt128 before; + while (!carryBits.IsZero) + { + before = vLow; + vLow.Add(carryBits); + carryBits = vLow.CompareTo(before) < 0 ? GpuUInt128.One : GpuUInt128.Zero; + } + + // Final correction if v >= modulus (mask) + if (vLow.High > maskHigh || (vLow.High == maskHigh && vLow.Low >= maskLow)) + { + ulong borrow = 0UL; + ulong newLow = vLow.Low - maskLow; + borrow = vLow.Low < maskLow ? 1UL : 0UL; + ulong newHigh = vLow.High - maskHigh - borrow; + vLow = new GpuUInt128(newHigh, newLow); + } + + return vLow; + } + + // Helper: 128-bit square -> 256-bit product as four 64-bit limbs (q3..q0) + private static void SquareFull128(GpuUInt128 value, out ulong q3, out ulong q2, out ulong q1, out ulong q0) + { + ulong L = value.Low; + ulong H = value.High; + Mul64Parts(L, L, out var hLL, out var lLL); // L^2 + Mul64Parts(H, H, out var hHH, out var lHH); // H^2 + Mul64Parts(L, H, out var hLH, out var lLH); // L*H + + // double LH + ulong dLH_low = lLH << 1; + ulong carry = (lLH >> 63) & 1UL; + ulong dLH_high = (hLH << 1) | carry; + + q0 = lLL; + ulong sum1 = hLL + dLH_low; + ulong c1 = sum1 < hLL ? 1UL : 0UL; + q1 = sum1; + + ulong sum2 = lHH + dLH_high; + ulong c2 = sum2 < lHH ? 1UL : 0UL; + sum2 += c1; + if (sum2 < c1) + { + c2++; + } + q2 = sum2; + q3 = hHH + c2; + } + + // 64x64 -> 128 multiply (high,low) + [MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] + private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low) + { + ulong a0 = (uint)a; + ulong a1 = a >> 32; + ulong b0 = (uint)b; + ulong b1 = b >> 32; + + ulong lo = a0 * b0; + ulong mid1 = a1 * b0; + ulong mid2 = a0 * b1; + ulong hi = a1 * b1; + + ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2; + low = (lo & 0xFFFFFFFFUL) | (carry << 32); + hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); + high = hi; + } + + [MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] + private static ulong ModPow64(ulong value, ulong exponent, ulong modulus) + { + ulong result = 1UL; + ulong baseValue = value % modulus; + ulong exp = exponent; + + while (exp != 0UL) + { + if ((exp & 1UL) != 0UL) + { + result = (ulong)(((UInt128)result * baseValue) % modulus); + } + + baseValue = (ulong)(((UInt128)baseValue * baseValue) % modulus); + exp >>= 1; + } + + return result; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/NttBarrettKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/NttBarrettKernels.cs new file mode 100644 index 0000000..1a4f700 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/NttBarrettKernels.cs @@ -0,0 +1,226 @@ +using System.Runtime.CompilerServices; +using ILGPU; + +namespace PerfectNumbers.Core.Gpu; + +internal static class NttBarrettKernels +{ + // Barrett-based scaling: v = v * scale mod n for 128-bit n + public static void ScaleBarrett128Kernel(Index1D index, ArrayView data, GpuUInt128 scale, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) + { + var v = data[index]; + // 128x128 -> 256 + Mul128Full(v.High, v.Low, scale.High, scale.Low, out var z3, out var z2, out var z1, out var z0); + // q approx + Mul3x2Top(z3, z2, z1, muHigh, muLow, out var qHigh, out var qLow); + // q*n + Mul128Full(qHigh, qLow, modHigh, modLow, out var qq3, out var qq2, out var qq1, out var qq0); + // r = z - q*n + Sub256(z3, z2, z1, z0, qq3, qq2, qq1, qq0, out var r3, out var r2, out var r1, out var r0); + ulong resHigh = r1; + ulong resLow = r0; + if (Geq128(resHigh, resLow, modHigh, modLow)) + { + Sub128(ref resHigh, ref resLow, modHigh, modLow); + if (Geq128(resHigh, resLow, modHigh, modLow)) + { + Sub128(ref resHigh, ref resLow, modHigh, modLow); + } + } + + data[index] = new GpuUInt128(resHigh, resLow); + } + + // Barrett-based squaring: v = v^2 mod n for 128-bit n + public static void SquareBarrett128Kernel(Index1D index, ArrayView data, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) + { + var v = data[index]; + Mul128Full(v.High, v.Low, v.High, v.Low, out var z3, out var z2, out var z1, out var z0); + Mul3x2Top(z3, z2, z1, muHigh, muLow, out var qHigh, out var qLow); + Mul128Full(qHigh, qLow, modHigh, modLow, out var qq3, out var qq2, out var qq1, out var qq0); + Sub256(z3, z2, z1, z0, qq3, qq2, qq1, qq0, out var r3, out var r2, out var r1, out var r0); + ulong resHigh = r1; + ulong resLow = r0; + if (Geq128(resHigh, resLow, modHigh, modLow)) + { + Sub128(ref resHigh, ref resLow, modHigh, modLow); + if (Geq128(resHigh, resLow, modHigh, modLow)) + { + Sub128(ref resHigh, ref resLow, modHigh, modLow); + } + } + data[index] = new GpuUInt128(resHigh, resLow); + } + + // Barrett reduction based stage for 128-bit moduli (no '%'). + public static void StageBarrett128Kernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddles, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) + { + int t = index.X; + // TODO: Use the bitmask-based remainder helper here as well to remove `%` from the butterfly stage and align with the + // optimized kernels highlighted in the GPU pow2mod benchmarks. + int j = t % half; + int block = t / half; + int k = block * len; + int i1 = k + j; + int i2 = i1 + half; + var u = data[i1]; + var v = data[i2]; + var w = twiddles[stageOffset + j]; + + // Compute v * w (128x128) -> 256-bit z3..z0 + Mul128Full(v.High, v.Low, w.High, w.Low, out var z3, out var z2, out var z1, out var z0); + + // t = floor(z / b) where b=2^64 => take (z3,z2,z1) + // Compute q = floor((t * mu) / b^3). Only need top two limbs of 5-limb product. + Mul3x2Top(z3, z2, z1, muHigh, muLow, out var qHigh, out var qLow); + + // q * n (128x128) -> 256-bit qq3..qq0 + Mul128Full(qHigh, qLow, modHigh, modLow, out var qq3, out var qq2, out var qq1, out var qq0); + + // r = z - q*n (256-bit) + Sub256(z3, z2, z1, z0, qq3, qq2, qq1, qq0, out var r3, out var r2, out var r1, out var r0); + + // Reduce r to [0, n) by at most two subtractions + // Discard higher limbs (should be zero if r < 2n), keep low 128-bit + ulong resHigh = r1; // r1 is limb1 after subtraction (since lower two limbs are r1 (high) and r0 (low)) + ulong resLow = r0; + + // while (res >= n) res -= n; (at most twice) + if (Geq128(resHigh, resLow, modHigh, modLow)) + { + Sub128(ref resHigh, ref resLow, modHigh, modLow); + if (Geq128(resHigh, resLow, modHigh, modLow)) + { + Sub128(ref resHigh, ref resLow, modHigh, modLow); + } + } + + var sum = new GpuUInt128(u); + var mulred = new GpuUInt128(resHigh, resLow); + sum.AddMod(mulred, new GpuUInt128(modHigh, modLow)); + u.SubMod(mulred, new GpuUInt128(modHigh, modLow)); + data[i1] = sum; + data[i2] = u; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Mul128Full(ulong aHigh, ulong aLow, ulong bHigh, ulong bLow, out ulong p3, out ulong p2, out ulong p1, out ulong p0) + { + // Full 128x128 -> 256 multiply using 64-bit limbs + Mul64Parts(aLow, bLow, out var h0, out var l0); + Mul64Parts(aLow, bHigh, out var h1, out var l1); + Mul64Parts(aHigh, bLow, out var h2, out var l2); + Mul64Parts(aHigh, bHigh, out var h3, out var l3); + + p0 = l0; + ulong carry = 0UL; + ulong sum = h0; + sum += l1; if (sum < l1) carry++; + sum += l2; if (sum < l2) carry++; + p1 = sum; + + sum = h1; + sum += h2; ulong carry2 = sum < h2 ? 1UL : 0UL; + sum += l3; if (sum < l3) carry2++; + sum += carry; if (sum < carry) carry2++; + p2 = sum; + p3 = h3 + carry2; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Mul3x2Top(ulong t2, ulong t1, ulong t0, ulong muHigh, ulong muLow, out ulong top1, out ulong top0) + { + // Compute top two limbs of ( (t2,t1,t0) * (muHigh,muLow) ) after shifting right by 192 bits. + // We build full 5 limbs and return p4 (top1) and p3 (top0). + ulong p0 = 0, p1 = 0, p2 = 0, p3 = 0, p4 = 0; + + void Add128To(ref ulong lo, ref ulong hi, ulong addLo, ulong addHi) + { + ulong s = lo + addLo; + ulong c = s < lo ? 1UL : 0UL; + lo = s; + ulong s2 = hi + addHi + c; + hi = s2; + } + + // t0 * muLow -> contributes to p0,p1 + Mul64Parts(t0, muLow, out var h, out var l); + Add128To(ref p0, ref p1, l, h); + + // t0 * muHigh -> to p1,p2 + Mul64Parts(t0, muHigh, out h, out l); + Add128To(ref p1, ref p2, l, h); + + // t1 * muLow -> to p1,p2 + Mul64Parts(t1, muLow, out h, out l); + Add128To(ref p1, ref p2, l, h); + + // t1 * muHigh -> to p2,p3 + Mul64Parts(t1, muHigh, out h, out l); + Add128To(ref p2, ref p3, l, h); + + // t2 * muLow -> to p2,p3 + Mul64Parts(t2, muLow, out h, out l); + Add128To(ref p2, ref p3, l, h); + + // t2 * muHigh -> to p3,p4 + Mul64Parts(t2, muHigh, out h, out l); + Add128To(ref p3, ref p4, l, h); + + top1 = p4; top0 = p3; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Sub256(ulong a3, ulong a2, ulong a1, ulong a0, ulong b3, ulong b2, ulong b1, ulong b0, out ulong r3, out ulong r2, out ulong r1, out ulong r0) + { + ulong borrow = 0; + r0 = SubWithBorrow(a0, b0, ref borrow); + r1 = SubWithBorrow(a1, b1, ref borrow); + r2 = SubWithBorrow(a2, b2, ref borrow); + r3 = SubWithBorrow(a3, b3, ref borrow); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong SubWithBorrow(ulong a, ulong b, ref ulong borrow) + { + ulong res = a - b - borrow; + borrow = ((a ^ b) & (a ^ res) & 0x8000_0000_0000_0000UL) != 0 ? 1UL : (a < b + borrow ? 1UL : 0UL); + return res; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool Geq128(ulong aHigh, ulong aLow, ulong bHigh, ulong bLow) + { + return aHigh > bHigh || (aHigh == bHigh && aLow >= bLow); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Sub128(ref ulong aHigh, ref ulong aLow, ulong bHigh, ulong bLow) + { + ulong borrow = 0; + ulong low = aLow - bLow; + borrow = aLow < bLow ? 1UL : 0UL; + ulong high = aHigh - bHigh - borrow; + aHigh = high; aLow = low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low) + { + // 64x64 -> 128 using 32-bit limbs + ulong a0 = (uint)a; + ulong a1 = a >> 32; + ulong b0 = (uint)b; + ulong b1 = b >> 32; + + ulong lo = a0 * b0; + ulong mid1 = a1 * b0; + ulong mid2 = a0 * b1; + ulong hi = a1 * b1; + + ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2; + low = (lo & 0xFFFFFFFFUL) | (carry << 32); + hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); + high = hi; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/NttButterflyKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/NttButterflyKernels.cs new file mode 100644 index 0000000..7e60e7e --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/NttButterflyKernels.cs @@ -0,0 +1,30 @@ +using ILGPU; + +namespace PerfectNumbers.Core.Gpu; + +internal static class NttButterflyKernels +{ + public static void StageKernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddles, GpuUInt128 modulus) + { + // TODO(NTT-OPT): Consider shared-memory tiling to improve locality: + // load a block of size `len` into shared memory, do butterflies, write back. + // Requires explicit grouped kernels and chosen group size. + int t = index.X; + // TODO: Replace this `%` with a bitmask when `half` is a power of two so the stage index math matches the faster + // residue strategy from the GPU residue benchmarks. + int j = t % half; + int block = t / half; + int k = block * len; + int i1 = k + j; + int i2 = i1 + half; + var u = data[i1]; + var v = data[i2]; + var w = twiddles[stageOffset + j]; + v.MulMod(w, modulus); + var sum = new GpuUInt128(u); + sum.AddMod(v, modulus); + u.SubMod(v, modulus); + data[i1] = sum; + data[i2] = u; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/NttMontgomeryKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/NttMontgomeryKernels.cs new file mode 100644 index 0000000..c58108d --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/NttMontgomeryKernels.cs @@ -0,0 +1,109 @@ +using System.Runtime.CompilerServices; +using ILGPU; + +namespace PerfectNumbers.Core.Gpu; + +internal static class NttMontgomeryKernels +{ + // TODO(MOD-OPT): Stage kernel using 64-bit Montgomery multiplication. Requires data and twiddles in Montgomery domain. + public static void StageMontKernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddlesMont, ulong modulus, ulong nPrime) + { + int t = index.X; + // TODO: Apply the bitmask remainder helper to this stage too so every butterfly path drops the slower `%` operation. + int j = t % half; + int block = t / half; + int k = block * len; + int i1 = k + j; + int i2 = i1 + half; + ulong u = data[i1].Low; + ulong v = data[i2].Low; + ulong wR = twiddlesMont[stageOffset + j].Low; + v = MontMul64(v, wR, modulus, nPrime); + ulong sum = u + v; + if (sum >= modulus) + { + sum -= modulus; + } + ulong diff = u >= v ? u - v : unchecked(u + modulus - v); + data[i1] = new GpuUInt128(0UL, sum); + data[i2] = new GpuUInt128(0UL, diff); + } + + public static void ToMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime, ulong r2) + { + ulong a = data[index].Low; + ulong v = MontMul64(a, r2, modulus, nPrime); + data[index] = new GpuUInt128(0UL, v); + } + + public static void FromMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime) + { + ulong aR = data[index].Low; + ulong v = MontMul64(aR, 1UL, modulus, nPrime); + data[index] = new GpuUInt128(0UL, v); + } + + public static void SquareMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime) + { + ulong aR = data[index].Low; + ulong v = MontMul64(aR, aR, modulus, nPrime); + data[index] = new GpuUInt128(0UL, v); + } + + public static void ScaleMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime, ulong scaleMont) + { + ulong aR = data[index].Low; + ulong v = MontMul64(aR, scaleMont, modulus, nPrime); + data[index] = new GpuUInt128(0UL, v); + } + + // Computes Montgomery multiplication for 64-bit operands in Montgomery domain. + // Returns a*b*R^{-1} mod modulus, where R=2^64 and nPrime = -modulus^{-1} mod 2^64. + private static ulong MontMul64(ulong aR, ulong bR, ulong modulus, ulong nPrime) + { + // t = aR * bR (128-bit) + ulong tLow, tHigh; + Mul64Parts(aR, bR, out tHigh, out tLow); + // m = (tLow * nPrime) mod 2^64 -> just low 64 bits of the product + ulong mLow, mHigh; + Mul64Parts(tLow, nPrime, out mHigh, out mLow); + // u = (t + m*modulus) >> 64 + ulong mmLow, mmHigh; + Mul64Parts(mLow, modulus, out mmHigh, out mmLow); + // add tLow + mmLow, keep carry + ulong carry = 0UL; + ulong low = tLow + mmLow; + if (low < tLow) + { + carry = 1UL; + } + // high = tHigh + mmHigh + carry + ulong high = tHigh + mmHigh + carry; + ulong u = high; // this is (t + m*n) >> 64 + if (u >= modulus) + { + u -= modulus; + } + return u; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low) + { + // 64x64 -> 128 using 32-bit limbs + ulong a0 = (uint)a; + ulong a1 = a >> 32; + ulong b0 = (uint)b; + ulong b1 = b >> 32; + + ulong lo = a0 * b0; + ulong mid1 = a1 * b0; + ulong mid2 = a0 * b1; + ulong hi = a1 * b1; + + ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2; + low = (lo & 0xFFFFFFFFUL) | (carry << 32); + hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); + high = hi; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/NttPointwiseKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/NttPointwiseKernels.cs new file mode 100644 index 0000000..66a6e28 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/NttPointwiseKernels.cs @@ -0,0 +1,20 @@ +using ILGPU; + +namespace PerfectNumbers.Core.Gpu; + +internal static class NttPointwiseKernels +{ + public static void MulKernel(Index1D index, ArrayView a, ArrayView b, GpuUInt128 modulus) + { + var val = a[index]; + val.MulMod(b[index], modulus); + a[index] = val; + } + + public static void ScaleKernel(Index1D index, ArrayView data, GpuUInt128 scale, GpuUInt128 modulus) + { + var v = data[index]; + v.MulMod(scale, modulus); + data[index] = v; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/NttTransformKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/NttTransformKernels.cs new file mode 100644 index 0000000..2740246 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/NttTransformKernels.cs @@ -0,0 +1,69 @@ +using ILGPU; + +namespace PerfectNumbers.Core.Gpu; + +internal static class NttTransformKernels +{ + public static void ForwardKernel(Index1D index, ArrayView input, ArrayView output, int length, GpuUInt128 modulus, GpuUInt128 root) + { + // TODO(NTT-OPT): Replace this reference O(n^2) kernel with a stage-wise + // Cooley–Tukey butterfly kernel (O(n log n)). This version uses a per- + // element ModPow in the hot loop which is far too slow and causes long + // single-kernel runtimes. After twiddle precomputation, each butterfly + // should perform: (u, v) -> (u+v*w, u-v*w) with a single MulMod. + var sum = new GpuUInt128(0UL, 0UL); + for (int i = 0; i < length; i++) + { + var term = input[i]; + var power = new GpuUInt128(root); + power.ModPow((ulong)index.X * (ulong)i, modulus); + term.MulMod(power, modulus); + sum.AddMod(term, modulus); + } + + output[index] = sum; + } + + public static void InverseKernel(Index1D index, ArrayView input, ArrayView output, int length, GpuUInt128 modulus, GpuUInt128 rootInv, GpuUInt128 nInv) + { + // TODO(NTT-OPT): Replace this O(n^2) inverse with stage-wise butterflies + // using precomputed inverse twiddles. Apply final scaling by nInv after + // all stages, preferably in a small dedicated kernel. + var sum = new GpuUInt128(0UL, 0UL); + for (int i = 0; i < length; i++) + { + var term = input[i]; + var power = new GpuUInt128(rootInv); + power.ModPow((ulong)index.X * (ulong)i, modulus); + term.MulMod(power, modulus); + sum.AddMod(term, modulus); + } + + sum.MulMod(nInv, modulus); + output[index] = sum; + } + + public static void BitReverseKernel(Index1D index, ArrayView data, int bits) + { + int i = index.X; + int j = ReverseBits(i, bits); + if (i < j) + { + var tmp = data[i]; + data[i] = data[j]; + data[j] = tmp; + } + } + + internal static int ReverseBits(int value, int bits) + { + int result = 0; + for (int i = 0; i < bits; i++) + { + result = (result << 1) | (value & 1); + value >>= 1; + } + + return result; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/NumericProbeKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/NumericProbeKernels.cs new file mode 100644 index 0000000..370632a --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/NumericProbeKernels.cs @@ -0,0 +1,51 @@ +using System.Numerics; +using ILGPU; +using ILGPU.Runtime; +using PeterO.Numbers; +using Open.Numeric.Primes; + +namespace PerfectNumbers.Core.Gpu; + +internal static class NumericProbeKernels +{ + public static void OpenNumericPrimesIsPrimeKernel(Index1D index, ArrayView1D data) + { + ulong value = (ulong)(index + 2); + bool isPrime = Prime.Numbers.IsPrime(value); + data[index] = isPrime ? 1 : 0; + } + + public static void BigIntegerKernel(Index1D index, ArrayView1D data) + { + BigInteger value = new BigInteger(data[index]); + value += BigInteger.One; + value -= BigInteger.One; + value *= new BigInteger(2); + value /= new BigInteger(3); + value %= new BigInteger(5); + int comparison = value.CompareTo(BigInteger.Zero); + data[index] = comparison; + } + + public static void EIntegerKernel(Index1D index, ArrayView1D data) + { + EInteger value = EInteger.FromInt32(data[index]); + value = value + EInteger.One; + value = value - EInteger.One; + value = value * EInteger.FromInt32(2); + value = value.Divide(EInteger.FromInt32(3)); + value = value.Remainder(EInteger.FromInt32(5)); + int comparison = value.CompareTo(EInteger.Zero); + data[index] = comparison; + } + + public static void ERationalKernel(Index1D index, ArrayView1D data) + { + ERational value = ERational.FromInt32(data[index]); + ERational adjusted = (value + ERational.FromInt32(1)) - ERational.FromInt32(1); + ERational product = adjusted * ERational.FromInt32(2); + ERational quotient = product / ERational.FromInt32(3); + int comparison = quotient.CompareTo(ERational.FromInt32(0)); + data[index] = comparison; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/OrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/OrderKernels.cs new file mode 100644 index 0000000..e17bdd1 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/OrderKernels.cs @@ -0,0 +1,51 @@ +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class OrderKernels +{ + public static void OrderKernelScan( + Index1D index, + ulong exponent, + ulong divMul, + ArrayView qs, + ArrayView orders) + { + GpuUInt128 q = qs[index]; + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + GpuUInt128 phi = q - GpuUInt128.One; + if (phi.High != 0UL) + { + orders[index] = 0UL; + return; + } + + ulong phi64 = phi.Low; + // TODO: Once the ProcessEightBitWindows helper is available, switch this order kernel to that faster + // Pow2Mod variant so cycle checks inherit the same gains observed in GpuPow2ModBenchmarks. + GpuUInt128 pow = GpuUInt128.Pow2Mod(phi64, in readOnlyQ); + if (pow != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + ulong div = KernelMathHelpers.FastDiv64Gpu(phi64, exponent, divMul); + GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, in readOnlyQ) - GpuUInt128.One; + if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) + { + orders[index] = 0UL; + return; + } + + orders[index] = exponent; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs new file mode 100644 index 0000000..85f06a4 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs @@ -0,0 +1,182 @@ +using ILGPU; +using ILGPU.Algorithms; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class Pow2ModKernels +{ + public static void Pow2ModKernelScan( + Index1D index, + ulong exponent, + GpuUInt128 twoP, + GpuUInt128 kStart, + byte lastIsSeven, + ulong _, + ResidueAutomatonArgs ra, + ArrayView orders, + ArrayView1D smallCycles, + ArrayView1D smallPrimesLastOne, + ArrayView1D smallPrimesLastSeven, + ArrayView1D smallPrimesPow2LastOne, + ArrayView1D smallPrimesPow2LastSeven) + { + ulong idx = (ulong)index.X; + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + orders[index] = 0UL; + return; + } + + // TODO: kStart is modified after this. Is this expected? + kStart.Add(idx); + GpuUInt128 q = twoP; + q.Mul64(kStart); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + // cycle should always be initialized if we're within array limit in production code + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + orders[index] = 0UL; + return; + } + } + // TODO: Swap Pow2Minus1Mod for the eight-bit window helper once the scalar version switches; + // benchmarks show the windowed variant cuts large-divisor scans from ~51 µs to ~21 µs. + if (GpuUInt128.Pow2Minus1Mod(exponent, in readOnlyQ) != GpuUInt128.Zero) + { + orders[index] = 0UL; + return; + } + + ArrayView1D primes = lastIsSeven != 0 ? smallPrimesLastSeven : smallPrimesLastOne; + ArrayView1D primesPow2 = lastIsSeven != 0 ? smallPrimesPow2LastSeven : smallPrimesPow2LastOne; + int primesLen = (int)primes.Length; + for (int i = 0; i < primesLen; i++) + { + ulong square = primesPow2[i]; + if (new GpuUInt128(0UL, square) > q) + { + break; + } + ulong prime = primes[i]; + if (KernelMathHelpers.Mod128By64(q, prime) == 0UL) + { + orders[index] = 0UL; + return; + } + } + + orders[index] = exponent; + } + + public static void Pow2ModOrderKernelScan( + Index1D index, + ulong exponent, + GpuUInt128 twoP, + GpuUInt128 kStart, + byte lastIsSeven, + ulong _, + ResidueAutomatonArgs ra, + ArrayView found, + ArrayView1D smallCycles) + { + ulong idx = (ulong)index.X; + ulong idxMod3 = idx % 3UL; + ulong idxMod5 = idx % 5UL; + ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + return; + } + + GpuUInt128 k = kStart + (GpuUInt128)idx; + GpuUInt128 q = twoP; + q.Mul64(k); + q.Add(GpuUInt128.One); + ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly(); + // Small-cycles in-kernel early rejection from device table + if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) + { + ulong cycle = smallCycles[(int)q.Low]; + if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) + { + return; + } + } + // TODO: Replace this Pow2Mod check with the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; + // residue order scans will then benefit from the same 2× speedup the benchmarked windowed kernel delivered. + if (GpuUInt128.Pow2Mod(exponent, in readOnlyQ) != GpuUInt128.One) + { + return; + } + + Atomic.Or(ref found[0], 1); + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs new file mode 100644 index 0000000..806e828 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs @@ -0,0 +1,17 @@ +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class Pow2MontgomeryKernels +{ + public static void Pow2MontgomeryKernel( + Index1D index, + ArrayView1D exponents, + MontgomeryDivisorData divisor, + byte keepMontgomery, + ArrayView1D results) + { + results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0); + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs new file mode 100644 index 0000000..493aa1b --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs @@ -0,0 +1,1358 @@ +using System; +using System.Numerics; +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Algorithms.Random; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class PrimeOrderKernels +{ + internal static void PartialFactorKernel( + Index1D index, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + int slotCount, + ulong value, + uint limit, + ArrayView1D factorsOut, + ArrayView1D exponentsOut, + ArrayView1D countOut, + ArrayView1D remainingOut, + ArrayView1D fullyFactoredOut) + { + if (index != 0) + { + return; + } + + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + ulong remainingLocal = value; + int count = 0; + + for (int i = 0; i < primeCount && count < slotCount; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remainingLocal) + { + break; + } + + ulong primeValue = primeCandidate; + if (primeValue == 0UL || (remainingLocal % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); + + factorsOut[count] = primeValue; + exponentsOut[count] = exponent; + count++; + } + + countOut[0] = count; + remainingOut[0] = remainingLocal; + fullyFactoredOut[0] = remainingLocal == 1UL ? (byte)1 : (byte)0; + } + + internal enum PrimeOrderKernelStatus : byte + { + Fallback = 0, + Found = 1, + HeuristicUnresolved = 2, + PollardOverflow = 3, + FactoringFailure = 4, + } + + private readonly struct CandidateKey + { + public CandidateKey(int primary, long secondary, long tertiary) + { + Primary = primary; + Secondary = secondary; + Tertiary = tertiary; + } + + public int Primary { get; } + + public long Secondary { get; } + + public long Tertiary { get; } + + public int CompareTo(CandidateKey other) + { + if (Primary != other.Primary) + { + return Primary.CompareTo(other.Primary); + } + + if (Secondary != other.Secondary) + { + return Secondary.CompareTo(other.Secondary); + } + + return Tertiary.CompareTo(other.Tertiary); + } + } + + internal static void CalculateOrderKernel( + Index1D index, + ulong prime, + OrderKernelConfig config, + MontgomeryDivisorData divisor, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + OrderKernelBuffers buffers) + { + if (index != 0) + { + return; + } + + ArrayView1D phiFactors = buffers.PhiFactors; + ArrayView1D phiExponents = buffers.PhiExponents; + ArrayView1D workFactors = buffers.WorkFactors; + ArrayView1D workExponents = buffers.WorkExponents; + ArrayView1D candidates = buffers.Candidates; + ArrayView1D stackIndex = buffers.StackIndex; + ArrayView1D stackExponent = buffers.StackExponent; + ArrayView1D stackProduct = buffers.StackProduct; + ArrayView1D resultOut = buffers.Result; + ArrayView1D statusOut = buffers.Status; + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong previousOrder = config.PreviousOrder; + byte hasPreviousOrder = config.HasPreviousOrder; + int maxPowChecks = config.MaxPowChecks; + int mode = config.Mode; + + statusOut[0] = (byte)PrimeOrderKernelStatus.Fallback; + + if (prime <= 3UL) + { + ulong orderValue = prime == 3UL ? 2UL : 1UL; + resultOut[0] = orderValue; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + ulong phi = prime - 1UL; + + int phiFactorCount = FactorWithSmallPrimes(phi, limit, primes, squares, primeCount, phiFactors, phiExponents, out ulong phiRemaining); + if (phiRemaining != 1UL) + { + // Reuse stackProduct as the Pollard Rho stack so factoring stays within this kernel. + if (!TryFactorWithPollardKernel( + phiRemaining, + limit, + primes, + squares, + primeCount, + phiFactors, + phiExponents, + ref phiFactorCount, + stackProduct, + statusOut)) + { + if (statusOut[0] != (byte)PrimeOrderKernelStatus.PollardOverflow) + { + resultOut[0] = CalculateByDoublingKernel(prime); + } + + return; + } + } + + SortFactors(phiFactors, phiExponents, phiFactorCount); + + if (TrySpecialMaxKernel(phi, prime, phiFactors, phiFactorCount, divisor)) + { + resultOut[0] = phi; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + ulong candidateOrder = InitializeStartingOrderKernel(prime, phi, divisor); + candidateOrder = ExponentLoweringKernel(candidateOrder, prime, phiFactors, phiExponents, phiFactorCount, divisor); + + if (TryConfirmOrderKernel( + prime, + candidateOrder, + divisor, + limit, + primes, + squares, + primeCount, + workFactors, + workExponents, + stackProduct, + statusOut)) + { + resultOut[0] = candidateOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + bool isStrict = mode == 1; + if (isStrict) + { + ulong strictOrder = CalculateByDoublingKernel(prime); + resultOut[0] = strictOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + if (TryHeuristicFinishKernel( + prime, + candidateOrder, + previousOrder, + hasPreviousOrder, + divisor, + limit, + maxPowChecks, + primes, + squares, + primeCount, + workFactors, + workExponents, + candidates, + stackIndex, + stackExponent, + stackProduct, + statusOut, + out ulong confirmedOrder)) + { + resultOut[0] = confirmedOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.Found; + return; + } + + ulong fallbackOrder = CalculateByDoublingKernel(prime); + resultOut[0] = fallbackOrder; + statusOut[0] = (byte)PrimeOrderKernelStatus.HeuristicUnresolved; + } + + private static int FactorWithSmallPrimes( + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + out ulong remaining) + { + remaining = value; + int factorCount = 0; + long factorLength = factors.Length; + long exponentLength = exponents.Length; + int capacity = factorLength < exponentLength ? (int)factorLength : (int)exponentLength; + + for (int i = 0; i < primeCount && remaining > 1UL && factorCount < capacity; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > limit) + { + break; + } + + ulong square = squares[i]; + if (square != 0UL && square > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + factors[factorCount] = primeValue; + exponents[factorCount] = exponent; + factorCount++; + } + + for (int i = factorCount; i < factors.Length; i++) + { + factors[i] = 0UL; + } + + for (int i = factorCount; i < exponents.Length; i++) + { + exponents[i] = 0; + } + + return factorCount; + } + + private static void SortFactors(ArrayView1D factors, ArrayView1D exponents, int count) + { + for (int i = 1; i < count; i++) + { + ulong factor = factors[i]; + int exponent = exponents[i]; + int j = i - 1; + + while (j >= 0 && factors[j] > factor) + { + factors[j + 1] = factors[j]; + exponents[j + 1] = exponents[j]; + j--; + } + + factors[j + 1] = factor; + exponents[j + 1] = exponent; + } + } + + private static bool TrySpecialMaxKernel( + ulong phi, + ulong prime, + ArrayView1D factors, + int factorCount, + in MontgomeryDivisorData divisor) + { + for (int i = 0; i < factorCount; i++) + { + ulong factor = factors[i]; + if (factor <= 1UL) + { + continue; + } + + ulong reduced = phi / factor; + if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL) + { + return false; + } + } + + return true; + } + + private static ulong InitializeStartingOrderKernel(ulong prime, ulong phi, in MontgomeryDivisorData divisor) + { + ulong order = phi; + ulong residue = prime & 7UL; + if (residue == 1UL || residue == 7UL) + { + ulong half = phi >> 1; + if (Pow2EqualsOneKernel(half, divisor)) + { + order = half; + } + } + + return order; + } + + private static ulong ExponentLoweringKernel( + ulong order, + ulong prime, + ArrayView1D factors, + ArrayView1D exponents, + int factorCount, + in MontgomeryDivisorData divisor) + { + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factors[i]; + int exponent = exponents[i]; + if (primeFactor <= 1UL) + { + continue; + } + + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong reduced = order / primeFactor; + if (Pow2EqualsOneKernel(reduced, divisor)) + { + order = reduced; + continue; + } + + break; + } + } + + return order; + } + + private static bool TryConfirmOrderKernel( + ulong prime, + ulong order, + in MontgomeryDivisorData divisor, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ArrayView1D compositeStack, + ArrayView1D statusOut) + { + if (order == 0UL) + { + return false; + } + + if (!Pow2EqualsOneKernel(order, divisor)) + { + return false; + } + + int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); + if (remaining != 1UL) + { + if (!TryFactorWithPollardKernel( + remaining, + limit, + primes, + squares, + primeCount, + factors, + exponents, + ref factorCount, + compositeStack, + statusOut)) + { + return false; + } + } + + SortFactors(factors, exponents, factorCount); + + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factors[i]; + int exponent = exponents[i]; + if (primeFactor <= 1UL) + { + continue; + } + + ulong reduced = order; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (Pow2EqualsOneKernel(reduced, divisor)) + { + return false; + } + } + } + + return true; + } + + private static bool TryHeuristicFinishKernel( + ulong prime, + ulong order, + ulong previousOrder, + byte hasPreviousOrder, + in MontgomeryDivisorData divisor, + uint limit, + int maxPowChecks, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D workFactors, + ArrayView1D workExponents, + ArrayView1D candidates, + ArrayView1D stackIndex, + ArrayView1D stackExponent, + ArrayView1D stackProduct, + ArrayView1D statusOut, + out ulong confirmedOrder) + { + confirmedOrder = 0UL; + + if (order <= 1UL) + { + return false; + } + + int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, workFactors, workExponents, out ulong remaining); + if (remaining != 1UL) + { + // Reuse stackProduct as the Pollard Rho stack while factoring the order candidates. + if (!TryFactorWithPollardKernel( + remaining, + limit, + primes, + squares, + primeCount, + workFactors, + workExponents, + ref factorCount, + stackProduct, + statusOut)) + { + return false; + } + } + + SortFactors(workFactors, workExponents, factorCount); + + long candidateCapacity = candidates.Length; + int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit; + int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit); + if (candidateCount == 0) + { + return false; + } + + SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); + + int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks; + if (powBudget <= 0) + { + powBudget = candidateCount; + } + + int powUsed = 0; + + for (int i = 0; i < candidateCount && powUsed < powBudget; i++) + { + ulong candidate = candidates[i]; + if (candidate <= 1UL) + { + continue; + } + + if (powUsed >= powBudget) + { + break; + } + + powUsed++; + if (!Pow2EqualsOneKernel(candidate, divisor)) + { + continue; + } + + if (!TryConfirmCandidateKernel( + prime, + candidate, + divisor, + limit, + primes, + squares, + primeCount, + workFactors, + workExponents, + stackProduct, + statusOut, + ref powUsed, + powBudget)) + { + continue; + } + + confirmedOrder = candidate; + return true; + } + + return false; + } + + private static int BuildCandidatesKernel( + ulong order, + ArrayView1D factors, + ArrayView1D exponents, + int factorCount, + ArrayView1D candidates, + ArrayView1D stackIndex, + ArrayView1D stackExponent, + ArrayView1D stackProduct, + int limit) + { + long stackIndexLength = stackIndex.Length; + long stackExponentLength = stackExponent.Length; + long stackProductLength = stackProduct.Length; + if (factorCount == 0 || limit <= 0 || stackIndexLength == 0L || stackExponentLength == 0L || stackProductLength == 0L) + { + return 0; + } + + int stackCapacity = stackIndexLength < stackExponentLength ? (int)stackIndexLength : (int)stackExponentLength; + if (stackProductLength < stackCapacity) + { + stackCapacity = (int)stackProductLength; + } + + int candidateCount = 0; + int stackTop = 0; + + stackIndex[0] = 0; + stackExponent[0] = 0; + stackProduct[0] = 1UL; + stackTop = 1; + + while (stackTop > 0) + { + stackTop--; + int index = stackIndex[stackTop]; + int exponent = stackExponent[stackTop]; + ulong product = stackProduct[stackTop]; + + if (index >= factorCount) + { + if (product != 1UL && product != order && candidateCount < limit) + { + ulong candidate = order / product; + if (candidate > 1UL && candidate < order) + { + candidates[candidateCount] = candidate; + candidateCount++; + } + } + + continue; + } + + int maxExponent = exponents[index]; + if (exponent > maxExponent) + { + continue; + } + + if (stackTop >= stackCapacity) + { + return candidateCount; + } + + stackIndex[stackTop] = index + 1; + stackExponent[stackTop] = 0; + stackProduct[stackTop] = product; + stackTop++; + + if (exponent == maxExponent) + { + continue; + } + + ulong primeFactor = factors[index]; + if (primeFactor == 0UL || product > order / primeFactor) + { + continue; + } + + if (stackTop >= stackCapacity) + { + return candidateCount; + } + + stackIndex[stackTop] = index; + stackExponent[stackTop] = exponent + 1; + stackProduct[stackTop] = product * primeFactor; + stackTop++; + } + + return candidateCount; + } + + private static void SortCandidatesKernel( + ulong prime, + ulong previousOrder, + bool hasPrevious, + ArrayView1D candidates, + int count) + { + for (int i = 1; i < count; i++) + { + ulong value = candidates[i]; + CandidateKey key = BuildCandidateKey(value, prime, previousOrder, hasPrevious); + int j = i - 1; + + while (j >= 0) + { + CandidateKey other = BuildCandidateKey(candidates[j], prime, previousOrder, hasPrevious); + if (other.CompareTo(key) <= 0) + { + break; + } + + candidates[j + 1] = candidates[j]; + j--; + } + + candidates[j + 1] = value; + } + } + + private static CandidateKey BuildCandidateKey(ulong value, ulong prime, ulong previousOrder, bool hasPrevious) + { + int group = GetGroup(value, prime); + if (group == 0) + { + return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); + } + + ulong reference = hasPrevious ? previousOrder : 0UL; + bool isGe = !hasPrevious || value >= reference; + int previousGroup = hasPrevious ? GetGroup(reference, prime) : 1; + int primary = ComputePrimary(group, isGe, previousGroup); + long secondary; + long tertiary; + + if (group == 3) + { + secondary = -(long)value; + tertiary = -(long)value; + } + else + { + ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; + secondary = (long)distance; + tertiary = (long)value; + } + + return new CandidateKey(primary, secondary, tertiary); + } + + private static int GetGroup(ulong value, ulong prime) + { + ulong threshold1 = prime >> 3; + if (value <= threshold1) + { + return 1; + } + + ulong threshold2 = prime >> 2; + if (value <= threshold2) + { + return 2; + } + + ulong threshold3 = (prime * 3UL) >> 3; + if (value <= threshold3) + { + return 3; + } + + return 0; + } + + private static int ComputePrimary(int group, bool isGe, int previousGroup) + { + int groupOffset; + switch (group) + { + case 1: + groupOffset = 0; + break; + case 2: + groupOffset = 2; + break; + case 3: + groupOffset = 4; + break; + default: + groupOffset = 6; + break; + } + + if (group == previousGroup) + { + if (group == 3) + { + return groupOffset + (isGe ? 0 : 3); + } + + return groupOffset + (isGe ? 0 : 1); + } + + return groupOffset + (isGe ? 0 : 1); + } + + private static bool TryConfirmCandidateKernel( + ulong prime, + ulong candidate, + in MontgomeryDivisorData divisor, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ArrayView1D compositeStack, + ArrayView1D statusOut, + ref int powUsed, + int powBudget) + { + int factorCount = FactorWithSmallPrimes(candidate, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); + if (remaining != 1UL) + { + if (!TryFactorWithPollardKernel( + remaining, + limit, + primes, + squares, + primeCount, + factors, + exponents, + ref factorCount, + compositeStack, + statusOut)) + { + return false; + } + } + + SortFactors(factors, exponents, factorCount); + + for (int i = 0; i < factorCount; i++) + { + ulong primeFactor = factors[i]; + int exponent = exponents[i]; + if (primeFactor <= 1UL) + { + continue; + } + + ulong reduced = candidate; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (powBudget > 0 && powUsed >= powBudget) + { + return false; + } + + powUsed++; + if (Pow2EqualsOneKernel(reduced, divisor)) + { + return false; + } + } + } + + return true; + } + + private static bool TryFactorWithPollardKernel( + ulong initial, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ref int factorCount, + ArrayView1D compositeStack, + ArrayView1D statusOut) + { + if (initial <= 1UL) + { + return true; + } + + int stackCapacity = (int)compositeStack.Length; + if (stackCapacity <= 0) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + int stackTop = 0; + compositeStack[stackTop] = initial; + stackTop++; + + while (stackTop > 0) + { + stackTop--; + ulong composite = compositeStack[stackTop]; + if (composite <= 1UL) + { + continue; + } + + if (!PeelSmallPrimesKernel( + composite, + limit, + primes, + squares, + primeCount, + factors, + exponents, + ref factorCount, + out ulong remaining)) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + if (remaining == 1UL) + { + continue; + } + + ulong factor = PollardRhoKernel(remaining); + if (factor <= 1UL || factor == remaining) + { + if (!TryAppendFactorKernel(factors, exponents, ref factorCount, remaining, 1)) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + continue; + } + + ulong quotient = remaining / factor; + if (stackTop + 2 > stackCapacity) + { + statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; + return false; + } + + compositeStack[stackTop] = factor; + compositeStack[stackTop + 1] = quotient; + stackTop += 2; + } + + return true; + } + + private static bool PeelSmallPrimesKernel( + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D factors, + ArrayView1D exponents, + ref int factorCount, + out ulong remaining) + { + ulong remainingLocal = value; + + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > limit) + { + break; + } + + ulong square = squares[i]; + if (square != 0UL && square > remainingLocal) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remainingLocal % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remainingLocal /= primeValue; + exponent++; + } + while ((remainingLocal % primeValue) == 0UL); + + if (!TryAppendFactorKernel(factors, exponents, ref factorCount, primeValue, exponent)) + { + remaining = value; + return false; + } + } + + remaining = remainingLocal; + return true; + } + + private static bool TryAppendFactorKernel( + ArrayView1D factors, + ArrayView1D exponents, + ref int count, + ulong prime, + int exponent) + { + if (prime <= 1UL || exponent <= 0) + { + return true; + } + + for (int i = 0; i < count; i++) + { + if (factors[i] == prime) + { + exponents[i] += exponent; + return true; + } + } + + int capacity = (int)factors.Length; + if (count >= capacity) + { + return false; + } + + factors[count] = prime; + exponents[count] = exponent; + count++; + return true; + } + + private static ulong MulModKernel(ulong left, ulong right, ulong modulus) + { + GpuUInt128 product = new GpuUInt128(left); + return product.MulMod(right, modulus); + } + + private static ulong PollardRhoKernel(ulong value) + { + if ((value & 1UL) == 0UL) + { + return 2UL; + } + + ulong c = 1UL; + while (true) + { + ulong x = 2UL; + ulong y = 2UL; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomialKernel(x, c, value); + y = AdvancePolynomialKernel(y, c, value); + y = AdvancePolynomialKernel(y, c, value); + + ulong diff = x > y ? x - y : y - x; + d = BinaryGcdKernel(diff, value); + } + + if (d == value) + { + c++; + if (c == 0UL) + { + c = 1UL; + } + + continue; + } + + return d; + } + } + + private static ulong AdvancePolynomialKernel(ulong x, ulong c, ulong modulus) + { + ulong squared = MulModKernel(x, x, modulus); + GpuUInt128 accumulator = new GpuUInt128(squared); + accumulator.AddMod(c, modulus); + return accumulator.Low; + } + + private static ulong BinaryGcdKernel(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + ulong aLocal = a >> BitOperations.TrailingZeroCount(a); + ulong bLocal = b; + + while (true) + { + bLocal >>= BitOperations.TrailingZeroCount(bLocal); + if (aLocal > bLocal) + { + ulong temp = aLocal; + aLocal = bLocal; + bLocal = temp; + } + + bLocal -= aLocal; + if (bLocal == 0UL) + { + return aLocal << shift; + } + } + } + + private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) + { + return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL; + } + + private static ulong CalculateByDoublingKernel(ulong prime) + { + ulong order = 1UL; + ulong value = 2UL % prime; + + while (value != 1UL) + { + value <<= 1; + if (value >= prime) + { + value -= prime; + } + + order++; + } + + return order; + } + + internal static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) + { + ulong exponent = exponents[index]; + remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false); + } + + private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) + { + return Pow2ModKernelWideCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(Pow2ModKernelWide); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); + }); + } + + internal static void Pow2ModKernelWide(Index1D index, ArrayView1D exponents, GpuUInt128 modulus, ArrayView1D remainders) + { + GpuUInt128 exponent = exponents[index]; + remainders[index] = Pow2ModKernelCore(exponent, modulus); + } + + private static GpuUInt128[] InitializeOddPowersTable(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + GpuUInt128[] result = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount]; + result[0] = baseValue; + if (oddPowerCount == 1) + { + return result; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + for (int i = 1; i < oddPowerCount; i++) + { + current.MulMod(baseValue, modulus); + result[i] = current; + } + + return result; + } + + internal static GpuUInt128 Pow2ModKernelCore(GpuUInt128 exponent, GpuUInt128 modulus) + { + // This should never happen in production code. + // if (modulus == GpuUInt128.One) + // { + // return GpuUInt128.Zero; + // } + + // This should never happen in production code. + // if (exponent.IsZero) + // { + // return GpuUInt128.One; + // } + + GpuUInt128 baseValue = UInt128Numbers.TwoGpu; + + // This should never happen in production code - 2 should never be greater or equal modulus. + // if (baseValue.CompareTo(modulus) >= 0) + // { + // baseValue.Sub(modulus); + // } + + if (ShouldUseSingleBit(exponent)) + { + return Pow2MontgomeryModSingleBit(exponent, modulus, baseValue); + } + + int index = exponent.GetBitLength(); + int windowSize = GetWindowSize(index); + index--; + + int oddPowerCount = 1 << (windowSize - 1); + GpuUInt128[] oddPowers = InitializeOddPowersTable(baseValue, modulus, oddPowerCount); + GpuUInt128 result = GpuUInt128.One; + + int windowStart; + ulong windowValue; + while (index >= 0) + { + if (!IsBitSet(exponent, index)) + { + result.MulMod(result, modulus); + index--; + continue; + } + + windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (!IsBitSet(exponent, windowStart)) + { + windowStart++; + } + + // We're reusing oddPowerCount as windowBitCount here to lower registry pressure & avoid additional allocation + oddPowerCount = index - windowStart + 1; + index = windowStart - 1; + windowValue = ExtractWindowValue(exponent, windowStart, oddPowerCount); + + // We're reusing windowStart as square here to lower registry pressure & avoid additional allocation + for (windowStart = 0; windowStart < oddPowerCount; windowStart++) + { + result.MulMod(result, modulus); + } + + // We're reusing windowStart as tableIndex here to lower registry pressure & avoid additional allocation + windowStart = (int)((windowValue - 1UL) >> 1); + result.MulMod(oddPowers[windowStart], modulus); + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool ShouldUseSingleBit(GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; + + private static int GetWindowSize(int bitLength) + { + if (bitLength <= Pow2WindowSizeBits) + { + return Math.Max(bitLength, 1); + } + + if (bitLength <= 23) + { + return 4; + } + + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + return Pow2WindowSizeBits; + } + + private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) + { + GpuUInt128 result = GpuUInt128.One; + + while (!exponent.IsZero) + { + if ((exponent.Low & 1UL) != 0UL) + { + result.MulMod(baseValue, modulus); + } + + exponent.ShiftRight(1); + if (exponent.IsZero) + { + break; + } + + // Reusing baseValue to store the squared base for the next iteration. + baseValue.MulMod(baseValue, modulus); + } + + return result; + } + + private static bool IsBitSet(GpuUInt128 value, int bitIndex) + { + if (bitIndex >= 64) + { + return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; + } + + return ((value.Low >> bitIndex) & 1UL) != 0UL; + } + + private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) + { + if (windowStart != 0) + { + exponent.ShiftRight(windowStart); + ulong mask = (1UL << windowBitCount) - 1UL; + return exponent.Low & mask; + } + + ulong directMask = (1UL << windowBitCount) - 1UL; + return exponent.Low & directMask; + } + +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs new file mode 100644 index 0000000..f831614 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -0,0 +1,99 @@ +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Algorithms; + +namespace PerfectNumbers.Core.Gpu; + +internal static class PrimeTesterKernels +{ + // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise. + public static void SmallPrimeSieveKernel(Index1D index, ArrayView numbers, ArrayView smallPrimes, ArrayView results) + { + ulong n = numbers[index]; + if (n <= 3UL) + { + results[index] = (byte)(n >= 2UL ? 1 : 0); + return; + } + + if ((n & 1UL) == 0UL || (n > 5UL && (n % 5UL) == 0UL)) + { + // TODO: Replace the `% 5` branch with the GPU Mod5 helper once the sieve kernel can + // reuse the benchmarked bitmask reduction instead of modulo on every candidate. + results[index] = 0; + return; + } + + if (n.Mod10() == 1UL) + { + // Early reject special GCD heuristic with floor(log2 n) + ulong m = 63UL - (ulong)XMath.LeadingZeroCount(n); + if (BinaryGcdGpu(n, m) != 1UL) + { + results[index] = 0; + return; + } + } + + long len = smallPrimes.Length; + for (int i = 0; i < len; i++) + { + ulong p = smallPrimes[i]; + if (p * p > n) + { + break; + } + + if ((n % p) == 0UL) + { + // TODO: Route this modulo through the shared divisor-cycle cache once exposed to GPU kernels + // so batched sieves avoid per-prime `%` operations that profiling showed expensive. + results[index] = 0; + return; + } + } + + results[index] = 1; + } + + public static void SharesFactorKernel(Index1D index, ArrayView numbers, ArrayView results) + { + ulong n = numbers[index]; + ulong m = 63UL - (ulong)XMath.LeadingZeroCount(n); + ulong gcd = BinaryGcdGpu(n, m); + results[index] = gcd == 1UL ? (byte)0 : (byte)1; + } + + private static ulong BinaryGcdGpu(ulong u, ulong v) + { + // TODO: Replace this inline GPU binary GCD with the kernel extracted from + // GpuUInt128BinaryGcdBenchmarks via GpuKernelPool so device callers reuse the + // fully unrolled ladder instead of this branchy fallback. + if (u == 0UL) + { + return v; + } + + if (v == 0UL) + { + return u; + } + + int shift = XMath.TrailingZeroCount(u | v); + u >>= XMath.TrailingZeroCount(u); + + do + { + v >>= XMath.TrailingZeroCount(v); + if (u > v) + { + (u, v) = (v, u); + } + + v -= u; + } + while (v != 0UL); + + return u << shift; + } +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/SmallPrimeFactorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/SmallPrimeFactorKernels.cs new file mode 100644 index 0000000..0ec44eb --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/SmallPrimeFactorKernels.cs @@ -0,0 +1,66 @@ +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal static class SmallPrimeFactorKernels +{ + public static void SmallPrimeFactorKernelScan( + Index1D index, + ulong value, + uint limit, + ArrayView1D primes, + ArrayView1D squares, + int primeCount, + ArrayView1D primeSlots, + ArrayView1D exponentSlots, + ArrayView1D countSlot, + ArrayView1D remainingSlot) + { + if (index != 0) + { + return; + } + + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + ulong remaining = value; + int count = 0; + int capacity = (int)primeSlots.Length; + + for (int i = 0; i < primeCount && remaining > 1UL && count < capacity; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate == 0U || primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + if (primeSquare != 0UL && primeSquare > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + primeSlots[count] = primeValue; + exponentSlots[count] = exponent; + count++; + } + + countSlot[0] = count; + remainingSlot[0] = remaining; + } +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 909176e..5868edc 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -28,23 +28,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly ConcurrentBag _sessionPool = new(); private Action, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(CheckKernel)); + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(DivisorByDivisorKernels.CheckKernel)); private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel)); + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.ComputeMontgomeryExponentKernel)); private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) => - _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(ComputeRemainderDeltasKernel)); + _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>(AccumulateRemaindersKernel)); + 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>(EvaluateCandidateMaskKernel)); + acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); public int GpuBatchSize { @@ -535,20 +535,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); } - private static void CheckKernel(Index1D index, ArrayView divisors, ArrayView exponents, ArrayView hits) - { - MontgomeryDivisorData divisor = divisors[index]; - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - - ulong exponent = exponents[index]; - hits[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; - } - public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession { private readonly MersenneNumberDivisorByDivisorGpuTester _owner; @@ -636,137 +622,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _owner._sessionPool.Add(this); } } - private static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) - { - deltas[index] = (byte)(gaps[index] % modulus); - } - - private static void AccumulateRemaindersKernel( - Index1D index, - ArrayView deltas, - ArrayView remainders, - byte baseRemainder, - byte modulus) - { - int globalIndex = index; - 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) - { - int addend = 0; - if (localIndex >= offset) - { - addend = shared[localIndex - offset]; - } - - Group.Barrier(); - - if (localIndex >= offset) - { - int sum = shared[localIndex] + addend; - if (sum >= mod) - { - sum %= mod; - } - - shared[localIndex] = sum; - } - - Group.Barrier(); - - offset <<= 1; - } - - if (globalIndex == 0) - { - remainders[0] = baseRemainder; - return; - } - - int remainder = baseRemainder + shared[localIndex]; - if (remainder >= mod) - { - remainder %= mod; - } - - remainders[globalIndex] = (byte)remainder; - } - - private static void EvaluateCandidateMaskKernel( - Index1D index, - ArrayView remainder10, - ArrayView remainder8, - ArrayView remainder5, - ArrayView remainder3, - byte lastIsSevenFlag, - ArrayView mask) - { - int globalIndex = index; - int length = (int)mask.Length; - if (globalIndex >= length) - { - return; - } - - byte value10 = remainder10[globalIndex]; - bool accept10; - if (lastIsSevenFlag != 0) - { - accept10 = value10 == 3 || value10 == 7 || value10 == 9; - } - else - { - accept10 = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9; - } - - if (!accept10) - { - mask[globalIndex] = 0; - return; - } - - byte value8 = remainder8[globalIndex]; - if (value8 != 1 && value8 != 7) - { - mask[globalIndex] = 0; - return; - } - - if (remainder3[globalIndex] == 0 || remainder5[globalIndex] == 0) - { - mask[globalIndex] = 0; - return; - } - - mask[globalIndex] = 1; - } - - private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) - { - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } - - ulong exponent = exponents[index]; - - results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true); - } public ulong DivisorLimit { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index 4ae660c..db1d688 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -1,6 +1,5 @@ using System.Collections.Concurrent; using ILGPU; -using ILGPU.Algorithms; using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu; @@ -14,7 +13,7 @@ public sealed class MersenneNumberDivisorGpuTester private readonly ConcurrentDictionary> _resultBuffers = new(); private Action> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(Kernel)); + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(DivisorKernels.Kernel)); public static void BuildDivisorCandidates() { @@ -131,47 +130,5 @@ public sealed class MersenneNumberDivisorGpuTester return true; } - private static void Kernel(Index1D _, ulong exponent, ReadOnlyGpuUInt128 divisor, ArrayView result) - { - if (divisor.High == 0UL && divisor.Low <= 1UL) - { - result[0] = 0; - return; - } - - GpuUInt128 baseVal; - int x = divisor.High == 0UL - ? 63 - XMath.LeadingZeroCount(divisor.Low) - : 127 - XMath.LeadingZeroCount(divisor.High); - ReadOnlyGpuUInt128 readOnlyMod = divisor; - if (x == 0) - { - // TODO: Swap Pow2Mod for the ProcessEightBitWindows helper once Pow2Minus1Mod adopts it; - // GpuPow2ModBenchmarks showed the windowed kernel cutting large divisors from ~51 µs to ~21 µs. - baseVal = GpuUInt128.Pow2Mod(exponent, in readOnlyMod); - result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0; - return; - } - - ulong ux = (ulong)x; - if (x >= 64) - { - baseVal = new GpuUInt128(1UL << (x - 64), 0UL); - } - else - { - baseVal = new GpuUInt128(0UL, 1UL << x); - } - - GpuUInt128 pow = baseVal; - pow.ModPow(exponent / ux, in readOnlyMod); - // TODO: Replace this trailing Pow2Mod call with the ProcessEightBitWindows variant once the shared - // windowed helper lands so mixed-radix decompositions benefit from the same GPU speedup. - var part2 = GpuUInt128.Pow2Mod(exponent % ux, in readOnlyMod); - GpuUInt128 product = pow; - ReadOnlyGpuUInt128 part2ReadOnly = part2.AsReadOnly(); - product.MulMod(in part2ReadOnly, in readOnlyMod); - result[0] = product.High == 0UL && product.Low == 1UL ? (byte)1 : (byte)0; - } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs index 76a137e..efac9d3 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs @@ -1,6 +1,5 @@ using System.Buffers; using System.Collections.Concurrent; -using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; using UInt128 = System.UInt128; @@ -18,23 +17,23 @@ public class MersenneNumberLucasLehmerGpuTester private Action> GetKernel(Accelerator accelerator) { - return KernelCache.GetOrAdd(accelerator, accel => accel.LoadAutoGroupedStreamKernel>(Kernel)); + return KernelCache.GetOrAdd(accelerator, accel => accel.LoadAutoGroupedStreamKernel>(LucasLehmerKernels.Kernel)); } private Action> GetAddSmallKernel(Accelerator accelerator) => - AddSmallKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(AddSmallKernel)); + AddSmallKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(LucasLehmerKernels.AddSmallKernel)); private Action> GetSubSmallKernel(Accelerator accelerator) => - SubSmallKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(SubtractSmallKernel)); + SubSmallKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(LucasLehmerKernels.SubtractSmallKernel)); private Action> GetReduceKernel(Accelerator accelerator) => - ReduceKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(ReduceModMersenneKernel)); + ReduceKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(LucasLehmerKernels.ReduceModMersenneKernel)); private Action, ArrayView> GetIsZeroKernel(Accelerator accelerator) => - IsZeroKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(IsZeroKernel)); + IsZeroKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(LucasLehmerKernels.IsZeroKernel)); private Action, ArrayView, ArrayView> GetBatchKernel(Accelerator accelerator) => - BatchKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(KernelBatch)); + BatchKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(LucasLehmerKernels.KernelBatch)); // Configurable LL slice size to keep kernels short. Default 32. public int SliceSize = 32; @@ -466,6 +465,13 @@ public class MersenneNumberLucasLehmerGpuTester return true; } + private static ulong MulMod(ulong a, ulong b, ulong modulus) + { + // TODO: Swap this UInt128 `%` reduction for the GPU-compatible MulMod helper once it adopts the faster + // inline UInt128 path benchmarked in MulMod64Benchmarks so host/GPU parity avoids BigInteger-style fallbacks. + return (ulong)(((UInt128)a * b) % modulus); + } + private static ulong ModPow(ulong value, ulong exponent, ulong modulus) { ulong result = 1UL; @@ -485,453 +491,5 @@ public class MersenneNumberLucasLehmerGpuTester return result; } - private static ulong MulMod(ulong a, ulong b, ulong modulus) - { - // TODO: Swap this UInt128 `%` reduction for the GPU-compatible MulMod helper once it adopts the faster - // inline UInt128 path benchmarked in MulMod64Benchmarks so host/GPU parity avoids BigInteger-style fallbacks. - return (ulong)(((UInt128)a * b) % modulus); - } - private static void AddSmallKernel(Index1D index, ulong add, ArrayView value) - { - ulong carry = add; - for (int i = 0; i < value.Length && carry != 0UL; i++) - { - var limb = value[i]; - ulong low = limb.Low + carry; - ulong carryOut = low < carry ? 1UL : 0UL; - ulong high = limb.High + carryOut; - carry = high < limb.High ? 1UL : 0UL; - value[i] = new GpuUInt128(high, low); - } - } - - private static void SubtractSmallKernel(Index1D index, ulong subtract, ArrayView value) - { - ulong borrow = subtract; - for (int i = 0; i < value.Length && borrow != 0UL; i++) - { - var limb = value[i]; - ulong low = limb.Low; - ulong high = limb.High; - if (low >= borrow) - { - low -= borrow; - borrow = 0UL; - } - else - { - low = unchecked(low - borrow); - borrow = 1UL; - if (high != 0UL) - { - high--; - borrow = 0UL; - } - else - { - high = ulong.MaxValue; - } - } - - value[i] = new GpuUInt128(high, low); - } - } - - private static void ReduceModMersenneKernel(Index1D index, ulong exponent, ArrayView value) - { - int limbCount = (int)((exponent + 127UL) / 128UL); - for (int i = limbCount; i < value.Length; i++) - { - var limb = value[i]; - if (limb.IsZero) - { - continue; - } - - int target = i - limbCount; - var sum = value[target]; - var original = sum; - sum.Add(limb); - value[target] = sum; - value[i] = new GpuUInt128(0UL, 0UL); - if (sum.CompareTo(original) < 0) - { - int j = target + 1; - while (j < limbCount) - { - var next = value[j]; - var nextOriginal = next; - next.Add(1UL); - value[j] = next; - if (next.CompareTo(nextOriginal) >= 0) - { - break; - } - - j++; - } - } - } - - int topBits = (int)exponent.Mod128(); - if (topBits == 0) - { - topBits = 128; - } - - int topIndex = limbCount - 1; - ulong topLowMask = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL; - ulong topHighMask = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL; - var top = value[topIndex]; - var carryBits = top >> topBits; - top = new(top.High & topHighMask, top.Low & topLowMask); - value[topIndex] = top; - - if (!carryBits.IsZero) - { - int j = 0; - while (j < limbCount && !carryBits.IsZero) - { - var cur = value[j]; - var curOriginal = cur; - cur.Add(carryBits); - value[j] = cur; - carryBits = cur.CompareTo(curOriginal) < 0 ? new GpuUInt128(0UL, 1UL) : new GpuUInt128(0UL, 0UL); - j++; - } - } - - bool geq = true; - for (int i = topIndex; i >= 0; i--) - { - var limb = value[i]; - ulong modHigh = i == topIndex ? topHighMask : ulong.MaxValue; - ulong modLow = i == topIndex ? topLowMask : ulong.MaxValue; - if (limb.High > modHigh || (limb.High == modHigh && limb.Low > modLow)) - { - geq = true; - break; - } - - if (limb.High < modHigh || (limb.High == modHigh && limb.Low < modLow)) - { - geq = false; - break; - } - } - - if (geq) - { - ulong borrow = 0UL; - for (int i = 0; i < limbCount - 1; i++) - { - var limb = value[i]; - ulong newLow = limb.Low - ulong.MaxValue - borrow; - borrow = limb.Low < ulong.MaxValue + borrow ? 1UL : 0UL; - ulong newHigh = limb.High - ulong.MaxValue - borrow; - borrow = limb.High < ulong.MaxValue + borrow ? 1UL : 0UL; - value[i] = new GpuUInt128(newHigh, newLow); - } - - var topLimb = value[topIndex]; - ulong newTopLow = topLimb.Low - topLowMask - borrow; - borrow = topLimb.Low < topLowMask + borrow ? 1UL : 0UL; - ulong newTopHigh = topLimb.High - topHighMask - borrow; - value[topIndex] = new GpuUInt128(newTopHigh, newTopLow); - } - - for (int i = limbCount; i < value.Length; i++) - { - value[i] = new GpuUInt128(0UL, 0UL); - } - } - - private static void IsZeroKernel(Index1D index, ArrayView value, ArrayView result) - { - byte isZero = 1; - for (int i = 0; i < value.Length; i++) - { - var limb = value[i]; - if (limb.High != 0UL || limb.Low != 0UL) - { - isZero = 0; - break; - } - } - - result[0] = isZero; - } - - private static void ReduceModMersenne(Span value, ulong exponent) - { - int limbCount = (int)((exponent + 127UL) / 128UL); - for (int i = limbCount; i < value.Length; i++) - { - var limb = value[i]; - if (limb.IsZero) - { - continue; - } - - int target = i - limbCount; - UInt128 original = value[target]; - UInt128 sum = original + (UInt128)limb; - bool carry = sum < original; - value[target] = new GpuUInt128(sum); - value[i] = new GpuUInt128(0UL, 0UL); - if (carry) - { - int j = target + 1; - while (j < limbCount) - { - UInt128 next = value[j]; - UInt128 nextSum = next + 1UL; - value[j] = new GpuUInt128(nextSum); - if (nextSum != 0UL) - { - break; - } - - j++; - } - } - } - - int topBits = (int)exponent.Mod128(); - if (topBits == 0) - { - topBits = 128; - } - - int topIndex = limbCount - 1; - ulong topLowMask = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL; - ulong topHighMask = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL; - UInt128 mask = ((UInt128)topHighMask << 64) | topLowMask; - UInt128 top = value[topIndex]; - UInt128 carryBits = top >> topBits; - top &= mask; - value[topIndex] = new GpuUInt128(top); - if (carryBits != 0UL) - { - int j = 0; - UInt128 carry = carryBits; - while (j < limbCount && carry != 0UL) - { - UInt128 cur = value[j]; - UInt128 sum = cur + carry; - bool overflow = sum < cur; - value[j] = new GpuUInt128(sum); - carry = overflow ? 1UL : 0UL; - j++; - } - } - - bool geq = true; - for (int i = topIndex; i >= 0; i--) - { - UInt128 limb = value[i]; - UInt128 modLimb = i == topIndex ? mask : UInt128.MaxValue; - if (limb > modLimb) - { - geq = true; - break; - } - - if (limb < modLimb) - { - geq = false; - break; - } - } - - if (geq) - { - ulong borrow = 0UL; - for (int i = 0; i < limbCount - 1; i++) - { - ulong low = value[i].Low; - ulong high = value[i].High; - ulong newLow = low - ulong.MaxValue - borrow; - borrow = low < ulong.MaxValue + borrow ? 1UL : 0UL; - ulong newHigh = high - ulong.MaxValue - borrow; - borrow = high < ulong.MaxValue + borrow ? 1UL : 0UL; - value[i] = new GpuUInt128(newHigh, newLow); - } - - ulong topLow = value[topIndex].Low; - ulong topHigh = value[topIndex].High; - ulong newTopLow = topLow - topLowMask - borrow; - borrow = topLow < topLowMask + borrow ? 1UL : 0UL; - ulong newTopHigh = topHigh - topHighMask - borrow; - value[topIndex] = new GpuUInt128(newTopHigh, newTopLow); - } - - for (int i = limbCount; i < value.Length; i++) - { - value[i] = new GpuUInt128(0UL, 0UL); - } - } - - private static void KernelBatch(Index1D index, ArrayView exponents, ArrayView moduli, ArrayView states) - { - int idx = index.X; - ulong exponent = exponents[idx]; - GpuUInt128 modulus = moduli[idx]; - var s = new GpuUInt128(0UL, 4UL); - ulong limit = exponent - 2UL; - for (ulong i = 0UL; i < limit; i++) - { - s = SquareModMersenne128(s, exponent); - s.SubMod(2UL, modulus); - } - - states[idx] = s; - } - - private static void Kernel(Index1D index, ulong exponent, GpuUInt128 modulus, ArrayView state) - { - // Lucas–Lehmer iteration in the field GF(2^p-1). - // For p < 128, use a Mersenne-specific squaring + reduction to avoid - // the generic 128-bit long-division-like reduction. - var s = new GpuUInt128(0UL, 4UL); - ulong i = 0UL, limit = exponent - 2UL; - if (exponent < 128UL) - { - for (; i < limit; i++) - { - s = SquareModMersenne128(s, exponent); - s.SubMod(2UL, modulus); - } - } - else - { - for (; i < limit; i++) - { - s.SquareMod(modulus); - s.SubMod(2UL, modulus); - } - } - - state[0] = s; - } - - // Squares a 128-bit value and reduces modulo M_p = 2^p - 1 using - // Mersenne folding. Valid for 0 < p < 128 and input s in [0, M_p]. - private static GpuUInt128 SquareModMersenne128(GpuUInt128 s, ulong p) - { - // Full square: (H*2^64 + L)^2 = H^2*2^128 + 2HL*2^64 + L^2 - SquareFull128(s, out var q3, out var q2, out var q1, out var q0); - - // v = low128 + high128 (fold by 128 bits once) - var vHigh = new GpuUInt128(q3, q2); - var vLow = new GpuUInt128(q1, q0); - vLow.Add(vHigh); - - // Mask to keep only p low bits and fold the carry bits (v >> p) - int topBits = (int)p.Mod128(); - if (topBits == 0) - { - topBits = 128; - } - - ulong maskLow = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL; - ulong maskHigh = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL; - - var carryBits = vLow >> topBits; - vLow = new(vLow.High & maskHigh, vLow.Low & maskLow); - - // Propagate carryBits back into the masked value (single-limb fold) - GpuUInt128 before; - while (!carryBits.IsZero) - { - before = vLow; - vLow.Add(carryBits); - carryBits = vLow.CompareTo(before) < 0 ? GpuUInt128.One : GpuUInt128.Zero; - } - - // Final correction if v >= modulus (mask) - if (vLow.High > maskHigh || (vLow.High == maskHigh && vLow.Low >= maskLow)) - { - ulong borrow = 0UL; - ulong newLow = vLow.Low - maskLow; - borrow = vLow.Low < maskLow ? 1UL : 0UL; - ulong newHigh = vLow.High - maskHigh - borrow; - vLow = new GpuUInt128(newHigh, newLow); - } - - return vLow; - } - - // Helper: 128-bit square -> 256-bit product as four 64-bit limbs (q3..q0) - private static void SquareFull128(GpuUInt128 value, out ulong q3, out ulong q2, out ulong q1, out ulong q0) - { - ulong L = value.Low; - ulong H = value.High; - Mul64Parts(L, L, out var hLL, out var lLL); // L^2 - Mul64Parts(H, H, out var hHH, out var lHH); // H^2 - Mul64Parts(L, H, out var hLH, out var lLH); // L*H - - // double LH - ulong dLH_low = lLH << 1; - ulong carry = (lLH >> 63) & 1UL; - ulong dLH_high = (hLH << 1) | carry; - - q0 = lLL; - ulong sum1 = hLL + dLH_low; - ulong c1 = sum1 < hLL ? 1UL : 0UL; - q1 = sum1; - - ulong sum2 = lHH + dLH_high; - ulong c2 = sum2 < lHH ? 1UL : 0UL; - sum2 += c1; - if (sum2 < c1) - { - c2++; - } - q2 = sum2; - q3 = hHH + c2; - } - - // 64x64 -> 128 multiply (high,low) - [MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] - private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low) - { - ulong a0 = (uint)a; - ulong a1 = a >> 32; - ulong b0 = (uint)b; - ulong b1 = b >> 32; - - ulong lo = a0 * b0; - ulong mid1 = a1 * b0; - ulong mid2 = a0 * b1; - ulong hi = a1 * b1; - - ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2; - low = (lo & 0xFFFFFFFFUL) | (carry << 32); - hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); - high = hi; - } - - [MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] - private static ulong ModPow64(ulong value, ulong exponent, ulong modulus) - { - ulong result = 1UL; - ulong baseValue = value % modulus; - ulong exp = exponent; - - while (exp != 0UL) - { - if ((exp & 1UL) != 0UL) - { - result = (ulong)(((UInt128)result * baseValue) % modulus); - } - - baseValue = (ulong)(((UInt128)baseValue * baseValue) % modulus); - exp >>= 1; - } - - return result; - } } - diff --git a/PerfectNumbers.Core/Gpu/NttGpuMath.cs b/PerfectNumbers.Core/Gpu/NttGpuMath.cs index 821a6b0..4f605ee 100644 --- a/PerfectNumbers.Core/Gpu/NttGpuMath.cs +++ b/PerfectNumbers.Core/Gpu/NttGpuMath.cs @@ -228,6 +228,26 @@ public static class NttGpuMath return 0UL - inv; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low) + { + // 64x64 -> 128 using 32-bit limbs + ulong a0 = (uint)a; + ulong a1 = a >> 32; + ulong b0 = (uint)b; + ulong b1 = b >> 32; + + ulong lo = a0 * b0; + ulong mid1 = a1 * b0; + ulong mid2 = a0 * b1; + ulong hi = a1 * b1; + + ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2; + low = (lo & 0xFFFFFFFFUL) | (carry << 32); + hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); + high = hi; + } + private readonly record struct SquareCacheKey(int Length, GpuUInt128 Modulus, GpuUInt128 PrimitiveRoot); private static readonly ConcurrentDictionary> SquareCache = new(); @@ -357,43 +377,43 @@ public static class NttGpuMath } private static Action, int> GetBitReverseKernel(Accelerator accelerator) => - BitReverseKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int>(BitReverseKernel)); + BitReverseKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int>(NttTransformKernels.BitReverseKernel)); private static Action, ArrayView, GpuUInt128> GetMulKernel(Accelerator accelerator) => - MulKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, GpuUInt128>(MulKernel)); + MulKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, GpuUInt128>(NttPointwiseKernels.MulKernel)); private static Action, int, int, int, ArrayView, GpuUInt128> GetStageKernel(Accelerator accelerator) => - StageKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int, int, int, ArrayView, GpuUInt128>(StageKernel)); + StageKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int, int, int, ArrayView, GpuUInt128>(NttButterflyKernels.StageKernel)); private static Action, GpuUInt128, GpuUInt128> GetScaleKernel(Accelerator accelerator) => - ScaleKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, GpuUInt128, GpuUInt128>(ScaleKernel)); + ScaleKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, GpuUInt128, GpuUInt128>(NttPointwiseKernels.ScaleKernel)); private static Action, int, int, int, ArrayView, ulong, ulong> GetStageMontKernel(Accelerator accelerator) => - StageMontKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int, int, int, ArrayView, ulong, ulong>(StageMontKernel)); + StageMontKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int, int, int, ArrayView, ulong, ulong>(NttMontgomeryKernels.StageMontKernel)); private static Action, int, int, int, ArrayView, ulong, ulong, ulong, ulong> GetStageBarrett128Kernel(Accelerator accelerator) => - StageBarrett128KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int, int, int, ArrayView, ulong, ulong, ulong, ulong>(StageBarrett128Kernel)); + StageBarrett128KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, int, int, int, ArrayView, ulong, ulong, ulong, ulong>(NttBarrettKernels.StageBarrett128Kernel)); private static Action, GpuUInt128, ulong, ulong, ulong, ulong> GetScaleBarrett128Kernel(Accelerator accelerator) => - ScaleBarrett128KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, GpuUInt128, ulong, ulong, ulong, ulong>(ScaleBarrett128Kernel)); + ScaleBarrett128KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, GpuUInt128, ulong, ulong, ulong, ulong>(NttBarrettKernels.ScaleBarrett128Kernel)); private static Action, ulong, ulong, ulong, ulong> GetSquareBarrett128Kernel(Accelerator accelerator) => - SquareBarrett128KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong>(SquareBarrett128Kernel)); + SquareBarrett128KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong>(NttBarrettKernels.SquareBarrett128Kernel)); private static Action, ulong, ulong, ulong> GetToMont64Kernel(Accelerator accelerator) => - ToMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong, ulong>(ToMont64Kernel)); + ToMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong, ulong>(NttMontgomeryKernels.ToMont64Kernel)); private static Action, ulong, ulong> GetFromMont64Kernel(Accelerator accelerator) => - FromMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong>(FromMont64Kernel)); + FromMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong>(NttMontgomeryKernels.FromMont64Kernel)); private static Action, ulong, ulong> GetSquareMont64Kernel(Accelerator accelerator) => - SquareMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong>(SquareMont64Kernel)); + SquareMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong>(NttMontgomeryKernels.SquareMont64Kernel)); private static Action, ulong, ulong, ulong> GetScaleMont64Kernel(Accelerator accelerator) => - ScaleMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong, ulong>(ScaleMont64Kernel)); + ScaleMont64KernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ulong, ulong, ulong>(NttMontgomeryKernels.ScaleMont64Kernel)); private static Action, ArrayView, int, GpuUInt128, GpuUInt128> GetForwardKernel(Accelerator accelerator) => - ForwardKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, int, GpuUInt128, GpuUInt128>(ForwardKernel)); + ForwardKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, int, GpuUInt128, GpuUInt128>(NttTransformKernels.ForwardKernel)); // TODO(NTT-OPT): Introduce stage-wise Cooley–Tukey kernels and a // per-(accelerator,length,modulus,root) twiddle cache. The plan: @@ -402,421 +422,19 @@ public static class NttGpuMath // - Provide accessors like GetStageKernel(...) and GetTwiddleBuffer(...). private static Action, ArrayView, int, GpuUInt128, GpuUInt128, GpuUInt128> GetInverseKernel(Accelerator accelerator) => - InverseKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, int, GpuUInt128, GpuUInt128, GpuUInt128>(InverseKernel)); + InverseKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, int, GpuUInt128, GpuUInt128, GpuUInt128>(NttTransformKernels.InverseKernel)); // TODO(NTT-OPT): Mirror the forward stage-wise design for the inverse // transform. Precompute inverse twiddles and a normalization factor // (n^-1 mod m). Launch one short kernel per stage to avoid TDR. - private static void MulKernel(Index1D index, ArrayView a, ArrayView b, GpuUInt128 modulus) - { - var val = a[index]; - val.MulMod(b[index], modulus); - a[index] = val; - } - - private static void StageKernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddles, GpuUInt128 modulus) - { - // TODO(NTT-OPT): Consider shared-memory tiling to improve locality: - // load a block of size `len` into shared memory, do butterflies, write back. - // Requires explicit grouped kernels and chosen group size. - int t = index.X; - // TODO: Replace this `%` with a bitmask when `half` is a power of two so the stage index math matches the faster - // residue strategy from the GPU residue benchmarks. - int j = t % half; - int block = t / half; - int k = block * len; - int i1 = k + j; - int i2 = i1 + half; - var u = data[i1]; - var v = data[i2]; - var w = twiddles[stageOffset + j]; - v.MulMod(w, modulus); - var sum = new GpuUInt128(u); - sum.AddMod(v, modulus); - u.SubMod(v, modulus); - data[i1] = sum; - data[i2] = u; - } - - private static void ScaleKernel(Index1D index, ArrayView data, GpuUInt128 scale, GpuUInt128 modulus) - { - var v = data[index]; - v.MulMod(scale, modulus); - data[index] = v; - } - - // Barrett-based scaling: v = v * scale mod n for 128-bit n - private static void ScaleBarrett128Kernel(Index1D index, ArrayView data, GpuUInt128 scale, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) - { - var v = data[index]; - // 128x128 -> 256 - Mul128Full(v.High, v.Low, scale.High, scale.Low, out var z3, out var z2, out var z1, out var z0); - // q approx - Mul3x2Top(z3, z2, z1, muHigh, muLow, out var qHigh, out var qLow); - // q*n - Mul128Full(qHigh, qLow, modHigh, modLow, out var qq3, out var qq2, out var qq1, out var qq0); - // r = z - q*n - Sub256(z3, z2, z1, z0, qq3, qq2, qq1, qq0, out var r3, out var r2, out var r1, out var r0); - ulong resHigh = r1; - ulong resLow = r0; - if (Geq128(resHigh, resLow, modHigh, modLow)) - { - Sub128(ref resHigh, ref resLow, modHigh, modLow); - if (Geq128(resHigh, resLow, modHigh, modLow)) - { - Sub128(ref resHigh, ref resLow, modHigh, modLow); - } - } - - data[index] = new GpuUInt128(resHigh, resLow); - } - - // Barrett-based squaring: v = v^2 mod n for 128-bit n - private static void SquareBarrett128Kernel(Index1D index, ArrayView data, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) - { - var v = data[index]; - Mul128Full(v.High, v.Low, v.High, v.Low, out var z3, out var z2, out var z1, out var z0); - Mul3x2Top(z3, z2, z1, muHigh, muLow, out var qHigh, out var qLow); - Mul128Full(qHigh, qLow, modHigh, modLow, out var qq3, out var qq2, out var qq1, out var qq0); - Sub256(z3, z2, z1, z0, qq3, qq2, qq1, qq0, out var r3, out var r2, out var r1, out var r0); - ulong resHigh = r1; - ulong resLow = r0; - if (Geq128(resHigh, resLow, modHigh, modLow)) - { - Sub128(ref resHigh, ref resLow, modHigh, modLow); - if (Geq128(resHigh, resLow, modHigh, modLow)) - { - Sub128(ref resHigh, ref resLow, modHigh, modLow); - } - } - data[index] = new GpuUInt128(resHigh, resLow); - } - - // Barrett reduction based stage for 128-bit moduli (no '%'). - private static void StageBarrett128Kernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddles, ulong modHigh, ulong modLow, ulong muHigh, ulong muLow) - { - int t = index.X; - // TODO: Use the bitmask-based remainder helper here as well to remove `%` from the butterfly stage and align with the - // optimized kernels highlighted in the GPU pow2mod benchmarks. - int j = t % half; - int block = t / half; - int k = block * len; - int i1 = k + j; - int i2 = i1 + half; - var u = data[i1]; - var v = data[i2]; - var w = twiddles[stageOffset + j]; - - // Compute v * w (128x128) -> 256-bit z3..z0 - Mul128Full(v.High, v.Low, w.High, w.Low, out var z3, out var z2, out var z1, out var z0); - - // t = floor(z / b) where b=2^64 => take (z3,z2,z1) - // Compute q = floor((t * mu) / b^3). Only need top two limbs of 5-limb product. - Mul3x2Top(z3, z2, z1, muHigh, muLow, out var qHigh, out var qLow); - - // q * n (128x128) -> 256-bit qq3..qq0 - Mul128Full(qHigh, qLow, modHigh, modLow, out var qq3, out var qq2, out var qq1, out var qq0); - - // r = z - q*n (256-bit) - Sub256(z3, z2, z1, z0, qq3, qq2, qq1, qq0, out var r3, out var r2, out var r1, out var r0); - - // Reduce r to [0, n) by at most two subtractions - // Discard higher limbs (should be zero if r < 2n), keep low 128-bit - ulong resHigh = r1; // r1 is limb1 after subtraction (since lower two limbs are r1 (high) and r0 (low)) - ulong resLow = r0; - - // while (res >= n) res -= n; (at most twice) - if (Geq128(resHigh, resLow, modHigh, modLow)) - { - Sub128(ref resHigh, ref resLow, modHigh, modLow); - if (Geq128(resHigh, resLow, modHigh, modLow)) - { - Sub128(ref resHigh, ref resLow, modHigh, modLow); - } - } - - var sum = new GpuUInt128(u); - var mulred = new GpuUInt128(resHigh, resLow); - sum.AddMod(mulred, new GpuUInt128(modHigh, modLow)); - u.SubMod(mulred, new GpuUInt128(modHigh, modLow)); - data[i1] = sum; - data[i2] = u; - } - - // TODO(MOD-OPT): Stage kernel using 64-bit Montgomery multiplication. Requires data and twiddles in Montgomery domain. - private static void StageMontKernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddlesMont, ulong modulus, ulong nPrime) - { - int t = index.X; - // TODO: Apply the bitmask remainder helper to this stage too so every butterfly path drops the slower `%` operation. - int j = t % half; - int block = t / half; - int k = block * len; - int i1 = k + j; - int i2 = i1 + half; - ulong u = data[i1].Low; - ulong v = data[i2].Low; - ulong wR = twiddlesMont[stageOffset + j].Low; - v = MontMul64(v, wR, modulus, nPrime); - ulong sum = u + v; - if (sum >= modulus) - { - sum -= modulus; - } - ulong diff = u >= v ? u - v : unchecked(u + modulus - v); - data[i1] = new GpuUInt128(0UL, sum); - data[i2] = new GpuUInt128(0UL, diff); - } - - // Computes Montgomery multiplication for 64-bit operands in Montgomery domain. - // Returns a*b*R^{-1} mod modulus, where R=2^64 and nPrime = -modulus^{-1} mod 2^64. - private static ulong MontMul64(ulong aR, ulong bR, ulong modulus, ulong nPrime) - { - // t = aR * bR (128-bit) - ulong tLow, tHigh; - Mul64Parts(aR, bR, out tHigh, out tLow); - // m = (tLow * nPrime) mod 2^64 -> just low 64 bits of the product - ulong mLow, mHigh; - Mul64Parts(tLow, nPrime, out mHigh, out mLow); - // u = (t + m*modulus) >> 64 - ulong mmLow, mmHigh; - Mul64Parts(mLow, modulus, out mmHigh, out mmLow); - // add tLow + mmLow, keep carry - ulong carry = 0UL; - ulong low = tLow + mmLow; - if (low < tLow) - { - carry = 1UL; - } - // high = tHigh + mmHigh + carry - ulong high = tHigh + mmHigh + carry; - ulong u = high; // this is (t + m*n) >> 64 - if (u >= modulus) - { - u -= modulus; - } - return u; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low) - { - // 64x64 -> 128 using 32-bit limbs - ulong a0 = (uint)a; - ulong a1 = a >> 32; - ulong b0 = (uint)b; - ulong b1 = b >> 32; - - ulong lo = a0 * b0; - ulong mid1 = a1 * b0; - ulong mid2 = a0 * b1; - ulong hi = a1 * b1; - - ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2; - low = (lo & 0xFFFFFFFFUL) | (carry << 32); - hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); - high = hi; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void Mul128Full(ulong aHigh, ulong aLow, ulong bHigh, ulong bLow, out ulong p3, out ulong p2, out ulong p1, out ulong p0) - { - // Full 128x128 -> 256 multiply using 64-bit limbs - Mul64Parts(aLow, bLow, out var h0, out var l0); - Mul64Parts(aLow, bHigh, out var h1, out var l1); - Mul64Parts(aHigh, bLow, out var h2, out var l2); - Mul64Parts(aHigh, bHigh, out var h3, out var l3); - - p0 = l0; - ulong carry = 0UL; - ulong sum = h0; - sum += l1; if (sum < l1) carry++; - sum += l2; if (sum < l2) carry++; - p1 = sum; - - sum = h1; - sum += h2; ulong carry2 = sum < h2 ? 1UL : 0UL; - sum += l3; if (sum < l3) carry2++; - sum += carry; if (sum < carry) carry2++; - p2 = sum; - p3 = h3 + carry2; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void Mul3x2Top(ulong t2, ulong t1, ulong t0, ulong muHigh, ulong muLow, out ulong top1, out ulong top0) - { - // Compute top two limbs of ( (t2,t1,t0) * (muHigh,muLow) ) after shifting right by 192 bits. - // We build full 5 limbs and return p4 (top1) and p3 (top0). - ulong p0 = 0, p1 = 0, p2 = 0, p3 = 0, p4 = 0; - - void Add128To(ref ulong lo, ref ulong hi, ulong addLo, ulong addHi) - { - ulong s = lo + addLo; - ulong c = s < lo ? 1UL : 0UL; - lo = s; - ulong s2 = hi + addHi + c; - hi = s2; - } - - // t0 * muLow -> contributes to p0,p1 - Mul64Parts(t0, muLow, out var h, out var l); - Add128To(ref p0, ref p1, l, h); - - // t0 * muHigh -> to p1,p2 - Mul64Parts(t0, muHigh, out h, out l); - Add128To(ref p1, ref p2, l, h); - - // t1 * muLow -> to p1,p2 - Mul64Parts(t1, muLow, out h, out l); - Add128To(ref p1, ref p2, l, h); - - // t1 * muHigh -> to p2,p3 - Mul64Parts(t1, muHigh, out h, out l); - Add128To(ref p2, ref p3, l, h); - - // t2 * muLow -> to p2,p3 - Mul64Parts(t2, muLow, out h, out l); - Add128To(ref p2, ref p3, l, h); - - // t2 * muHigh -> to p3,p4 - Mul64Parts(t2, muHigh, out h, out l); - Add128To(ref p3, ref p4, l, h); - - top1 = p4; top0 = p3; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void Sub256(ulong a3, ulong a2, ulong a1, ulong a0, ulong b3, ulong b2, ulong b1, ulong b0, out ulong r3, out ulong r2, out ulong r1, out ulong r0) - { - ulong borrow = 0; - r0 = SubWithBorrow(a0, b0, ref borrow); - r1 = SubWithBorrow(a1, b1, ref borrow); - r2 = SubWithBorrow(a2, b2, ref borrow); - r3 = SubWithBorrow(a3, b3, ref borrow); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong SubWithBorrow(ulong a, ulong b, ref ulong borrow) - { - ulong res = a - b - borrow; - borrow = ((a ^ b) & (a ^ res) & 0x8000_0000_0000_0000UL) != 0 ? 1UL : (a < b + borrow ? 1UL : 0UL); - return res; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool Geq128(ulong aHigh, ulong aLow, ulong bHigh, ulong bLow) - { - return aHigh > bHigh || (aHigh == bHigh && aLow >= bLow); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void Sub128(ref ulong aHigh, ref ulong aLow, ulong bHigh, ulong bLow) - { - ulong borrow = 0; - ulong low = aLow - bLow; - borrow = aLow < bLow ? 1UL : 0UL; - ulong high = aHigh - bHigh - borrow; - aHigh = high; aLow = low; - } - - private static void ToMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime, ulong r2) - { - ulong a = data[index].Low; - ulong v = MontMul64(a, r2, modulus, nPrime); - data[index] = new GpuUInt128(0UL, v); - } - - private static void FromMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime) - { - ulong aR = data[index].Low; - ulong v = MontMul64(aR, 1UL, modulus, nPrime); - data[index] = new GpuUInt128(0UL, v); - } - - private static void SquareMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime) - { - ulong aR = data[index].Low; - ulong v = MontMul64(aR, aR, modulus, nPrime); - data[index] = new GpuUInt128(0UL, v); - } - - private static void ScaleMont64Kernel(Index1D index, ArrayView data, ulong modulus, ulong nPrime, ulong scaleMont) - { - ulong aR = data[index].Low; - ulong v = MontMul64(aR, scaleMont, modulus, nPrime); - data[index] = new GpuUInt128(0UL, v); - } - - private static void ForwardKernel(Index1D index, ArrayView input, ArrayView output, int length, GpuUInt128 modulus, GpuUInt128 root) - { - // TODO(NTT-OPT): Replace this reference O(n^2) kernel with a stage-wise - // Cooley–Tukey butterfly kernel (O(n log n)). This version uses a per- - // element ModPow in the hot loop which is far too slow and causes long - // single-kernel runtimes. After twiddle precomputation, each butterfly - // should perform: (u, v) -> (u+v*w, u-v*w) with a single MulMod. - var sum = new GpuUInt128(0UL, 0UL); - for (int i = 0; i < length; i++) - { - var term = input[i]; - var power = new GpuUInt128(root); - power.ModPow((ulong)index.X * (ulong)i, modulus); - term.MulMod(power, modulus); - sum.AddMod(term, modulus); - } - - output[index] = sum; - } - - private static void InverseKernel(Index1D index, ArrayView input, ArrayView output, int length, GpuUInt128 modulus, GpuUInt128 rootInv, GpuUInt128 nInv) - { - // TODO(NTT-OPT): Replace this O(n^2) inverse with stage-wise butterflies - // using precomputed inverse twiddles. Apply final scaling by nInv after - // all stages, preferably in a small dedicated kernel. - var sum = new GpuUInt128(0UL, 0UL); - for (int i = 0; i < length; i++) - { - var term = input[i]; - var power = new GpuUInt128(rootInv); - power.ModPow((ulong)index.X * (ulong)i, modulus); - term.MulMod(power, modulus); - sum.AddMod(term, modulus); - } - - sum.MulMod(nInv, modulus); - output[index] = sum; - } - - private static void BitReverseKernel(Index1D index, ArrayView data, int bits) - { - int i = index.X; - int j = ReverseBits(i, bits); - if (i < j) - { - var tmp = data[i]; - data[i] = data[j]; - data[j] = tmp; - } - } - - private static int ReverseBits(int value, int bits) - { - int result = 0; - for (int i = 0; i < bits; i++) - { - result = (result << 1) | (value & 1); - value >>= 1; - } - - return result; - } - public static void BitReverse(Span values) { int n = values.Length; int bits = (int)Math.Log2(n); for (int i = 0; i < n; i++) { - int j = ReverseBits(i, bits); + int j = NttTransformKernels.ReverseBits(i, bits); if (i < j) { (values[j], values[i]) = (values[i], values[j]); diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 7ce5c57..9b77a3d 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -123,7 +123,7 @@ internal static partial class PrimeOrderGpuHeuristics { return PartialFactorKernelCache.GetOrAdd(accelerator, static accel => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(PartialFactorKernel); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(PrimeOrderKernels.PartialFactorKernel); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); }); @@ -266,66 +266,6 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static void PartialFactorKernel( - Index1D index, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - int slotCount, - ulong value, - uint limit, - ArrayView1D factorsOut, - ArrayView1D exponentsOut, - ArrayView1D countOut, - ArrayView1D remainingOut, - ArrayView1D fullyFactoredOut) - { - if (index != 0) - { - return; - } - - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - ulong remainingLocal = value; - int count = 0; - - for (int i = 0; i < primeCount && count < slotCount; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > effectiveLimit) - { - break; - } - - ulong primeSquare = squares[i]; - if (primeSquare != 0UL && primeSquare > remainingLocal) - { - break; - } - - ulong primeValue = primeCandidate; - if (primeValue == 0UL || (remainingLocal % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); - - factorsOut[count] = primeValue; - exponentsOut[count] = exponent; - count++; - } - - countOut[0] = count; - remainingOut[0] = remainingLocal; - fullyFactoredOut[0] = remainingLocal == 1UL ? (byte)1 : (byte)0; - } public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData divisorData) { Span exponents = stackalloc ulong[1]; @@ -480,15 +420,15 @@ internal static partial class PrimeOrderGpuHeuristics byte status = 0; statusBuffer.View.CopyToCPU(ref status, 1); - PrimeOrderKernelStatus kernelStatus = (PrimeOrderKernelStatus)status; - if (kernelStatus == PrimeOrderKernelStatus.Fallback) + PrimeOrderKernels.PrimeOrderKernelStatus kernelStatus = (PrimeOrderKernels.PrimeOrderKernelStatus)status; + if (kernelStatus == PrimeOrderKernels.PrimeOrderKernelStatus.Fallback) { DisposeResources(); lease.Dispose(); return false; } - if (kernelStatus == PrimeOrderKernelStatus.PollardOverflow) + if (kernelStatus == PrimeOrderKernels.PrimeOrderKernelStatus.PollardOverflow) { DisposeResources(); lease.Dispose(); @@ -497,7 +437,7 @@ internal static partial class PrimeOrderGpuHeuristics resultBuffer.View.CopyToCPU(ref order, 1); - if (kernelStatus == PrimeOrderKernelStatus.FactoringFailure) + if (kernelStatus == PrimeOrderKernels.PrimeOrderKernelStatus.FactoringFailure) { order = 0UL; } @@ -527,1094 +467,12 @@ internal static partial class PrimeOrderGpuHeuristics { return OrderKernelCache.GetOrAdd(accelerator, static accel => { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(CalculateOrderKernel); + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(PrimeOrderKernels.CalculateOrderKernel); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate(); }); } - private enum PrimeOrderKernelStatus : byte - { - Fallback = 0, - Found = 1, - HeuristicUnresolved = 2, - PollardOverflow = 3, - FactoringFailure = 4, - } - - private readonly struct CandidateKey - { - public CandidateKey(int primary, long secondary, long tertiary) - { - Primary = primary; - Secondary = secondary; - Tertiary = tertiary; - } - - public int Primary { get; } - - public long Secondary { get; } - - public long Tertiary { get; } - - public int CompareTo(CandidateKey other) - { - if (Primary != other.Primary) - { - return Primary.CompareTo(other.Primary); - } - - if (Secondary != other.Secondary) - { - return Secondary.CompareTo(other.Secondary); - } - - return Tertiary.CompareTo(other.Tertiary); - } - } - - private static void CalculateOrderKernel( - Index1D index, - ulong prime, - OrderKernelConfig config, - MontgomeryDivisorData divisor, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - OrderKernelBuffers buffers) - { - if (index != 0) - { - return; - } - - ArrayView1D phiFactors = buffers.PhiFactors; - ArrayView1D phiExponents = buffers.PhiExponents; - ArrayView1D workFactors = buffers.WorkFactors; - ArrayView1D workExponents = buffers.WorkExponents; - ArrayView1D candidates = buffers.Candidates; - ArrayView1D stackIndex = buffers.StackIndex; - ArrayView1D stackExponent = buffers.StackExponent; - ArrayView1D stackProduct = buffers.StackProduct; - ArrayView1D resultOut = buffers.Result; - ArrayView1D statusOut = buffers.Status; - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong previousOrder = config.PreviousOrder; - byte hasPreviousOrder = config.HasPreviousOrder; - int maxPowChecks = config.MaxPowChecks; - int mode = config.Mode; - - statusOut[0] = (byte)PrimeOrderKernelStatus.Fallback; - - if (prime <= 3UL) - { - ulong orderValue = prime == 3UL ? 2UL : 1UL; - resultOut[0] = orderValue; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - ulong phi = prime - 1UL; - - int phiFactorCount = FactorWithSmallPrimes(phi, limit, primes, squares, primeCount, phiFactors, phiExponents, out ulong phiRemaining); - if (phiRemaining != 1UL) - { - // Reuse stackProduct as the Pollard Rho stack so factoring stays within this kernel. - if (!TryFactorWithPollardKernel( - phiRemaining, - limit, - primes, - squares, - primeCount, - phiFactors, - phiExponents, - ref phiFactorCount, - stackProduct, - statusOut)) - { - if (statusOut[0] != (byte)PrimeOrderKernelStatus.PollardOverflow) - { - resultOut[0] = CalculateByDoublingKernel(prime); - } - - return; - } - } - - SortFactors(phiFactors, phiExponents, phiFactorCount); - - if (TrySpecialMaxKernel(phi, prime, phiFactors, phiFactorCount, divisor)) - { - resultOut[0] = phi; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - ulong candidateOrder = InitializeStartingOrderKernel(prime, phi, divisor); - candidateOrder = ExponentLoweringKernel(candidateOrder, prime, phiFactors, phiExponents, phiFactorCount, divisor); - - if (TryConfirmOrderKernel( - prime, - candidateOrder, - divisor, - limit, - primes, - squares, - primeCount, - workFactors, - workExponents, - stackProduct, - statusOut)) - { - resultOut[0] = candidateOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - bool isStrict = mode == 1; - if (isStrict) - { - ulong strictOrder = CalculateByDoublingKernel(prime); - resultOut[0] = strictOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - if (TryHeuristicFinishKernel( - prime, - candidateOrder, - previousOrder, - hasPreviousOrder, - divisor, - limit, - maxPowChecks, - primes, - squares, - primeCount, - workFactors, - workExponents, - candidates, - stackIndex, - stackExponent, - stackProduct, - statusOut, - out ulong confirmedOrder)) - { - resultOut[0] = confirmedOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.Found; - return; - } - - ulong fallbackOrder = CalculateByDoublingKernel(prime); - resultOut[0] = fallbackOrder; - statusOut[0] = (byte)PrimeOrderKernelStatus.HeuristicUnresolved; - } - - private static int FactorWithSmallPrimes( - ulong value, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - out ulong remaining) - { - remaining = value; - int factorCount = 0; - long factorLength = factors.Length; - long exponentLength = exponents.Length; - int capacity = factorLength < exponentLength ? (int)factorLength : (int)exponentLength; - - for (int i = 0; i < primeCount && remaining > 1UL && factorCount < capacity; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate == 0U || primeCandidate > limit) - { - break; - } - - ulong square = squares[i]; - if (square != 0UL && square > remaining) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == 0UL); - - factors[factorCount] = primeValue; - exponents[factorCount] = exponent; - factorCount++; - } - - for (int i = factorCount; i < factors.Length; i++) - { - factors[i] = 0UL; - } - - for (int i = factorCount; i < exponents.Length; i++) - { - exponents[i] = 0; - } - - return factorCount; - } - - private static void SortFactors(ArrayView1D factors, ArrayView1D exponents, int count) - { - for (int i = 1; i < count; i++) - { - ulong factor = factors[i]; - int exponent = exponents[i]; - int j = i - 1; - - while (j >= 0 && factors[j] > factor) - { - factors[j + 1] = factors[j]; - exponents[j + 1] = exponents[j]; - j--; - } - - factors[j + 1] = factor; - exponents[j + 1] = exponent; - } - } - - private static bool TrySpecialMaxKernel( - ulong phi, - ulong prime, - ArrayView1D factors, - int factorCount, - in MontgomeryDivisorData divisor) - { - for (int i = 0; i < factorCount; i++) - { - ulong factor = factors[i]; - if (factor <= 1UL) - { - continue; - } - - ulong reduced = phi / factor; - if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL) - { - return false; - } - } - - return true; - } - - private static ulong InitializeStartingOrderKernel(ulong prime, ulong phi, in MontgomeryDivisorData divisor) - { - ulong order = phi; - ulong residue = prime & 7UL; - if (residue == 1UL || residue == 7UL) - { - ulong half = phi >> 1; - if (Pow2EqualsOneKernel(half, divisor)) - { - order = half; - } - } - - return order; - } - - private static ulong ExponentLoweringKernel( - ulong order, - ulong prime, - ArrayView1D factors, - ArrayView1D exponents, - int factorCount, - in MontgomeryDivisorData divisor) - { - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factors[i]; - int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } - - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } - - ulong reduced = order / primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) - { - order = reduced; - continue; - } - - break; - } - } - - return order; - } - - private static bool TryConfirmOrderKernel( - ulong prime, - ulong order, - in MontgomeryDivisorData divisor, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ArrayView1D compositeStack, - ArrayView1D statusOut) - { - if (order == 0UL) - { - return false; - } - - if (!Pow2EqualsOneKernel(order, divisor)) - { - return false; - } - - int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); - if (remaining != 1UL) - { - if (!TryFactorWithPollardKernel( - remaining, - limit, - primes, - squares, - primeCount, - factors, - exponents, - ref factorCount, - compositeStack, - statusOut)) - { - return false; - } - } - - SortFactors(factors, exponents, factorCount); - - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factors[i]; - int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } - - ulong reduced = order; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) - { - return false; - } - } - } - - return true; - } - - private static bool TryHeuristicFinishKernel( - ulong prime, - ulong order, - ulong previousOrder, - byte hasPreviousOrder, - in MontgomeryDivisorData divisor, - uint limit, - int maxPowChecks, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D workFactors, - ArrayView1D workExponents, - ArrayView1D candidates, - ArrayView1D stackIndex, - ArrayView1D stackExponent, - ArrayView1D stackProduct, - ArrayView1D statusOut, - out ulong confirmedOrder) - { - confirmedOrder = 0UL; - - if (order <= 1UL) - { - return false; - } - - int factorCount = FactorWithSmallPrimes(order, limit, primes, squares, primeCount, workFactors, workExponents, out ulong remaining); - if (remaining != 1UL) - { - // Reuse stackProduct as the Pollard Rho stack while factoring the order candidates. - if (!TryFactorWithPollardKernel( - remaining, - limit, - primes, - squares, - primeCount, - workFactors, - workExponents, - ref factorCount, - stackProduct, - statusOut)) - { - return false; - } - } - - SortFactors(workFactors, workExponents, factorCount); - - long candidateCapacity = candidates.Length; - int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit; - int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit); - if (candidateCount == 0) - { - return false; - } - - SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); - - int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks; - if (powBudget <= 0) - { - powBudget = candidateCount; - } - - int powUsed = 0; - - for (int i = 0; i < candidateCount && powUsed < powBudget; i++) - { - ulong candidate = candidates[i]; - if (candidate <= 1UL) - { - continue; - } - - if (powUsed >= powBudget) - { - break; - } - - powUsed++; - if (!Pow2EqualsOneKernel(candidate, divisor)) - { - continue; - } - - if (!TryConfirmCandidateKernel( - prime, - candidate, - divisor, - limit, - primes, - squares, - primeCount, - workFactors, - workExponents, - stackProduct, - statusOut, - ref powUsed, - powBudget)) - { - continue; - } - - confirmedOrder = candidate; - return true; - } - - return false; - } - - private static int BuildCandidatesKernel( - ulong order, - ArrayView1D factors, - ArrayView1D exponents, - int factorCount, - ArrayView1D candidates, - ArrayView1D stackIndex, - ArrayView1D stackExponent, - ArrayView1D stackProduct, - int limit) - { - long stackIndexLength = stackIndex.Length; - long stackExponentLength = stackExponent.Length; - long stackProductLength = stackProduct.Length; - if (factorCount == 0 || limit <= 0 || stackIndexLength == 0L || stackExponentLength == 0L || stackProductLength == 0L) - { - return 0; - } - - int stackCapacity = stackIndexLength < stackExponentLength ? (int)stackIndexLength : (int)stackExponentLength; - if (stackProductLength < stackCapacity) - { - stackCapacity = (int)stackProductLength; - } - - int candidateCount = 0; - int stackTop = 0; - - stackIndex[0] = 0; - stackExponent[0] = 0; - stackProduct[0] = 1UL; - stackTop = 1; - - while (stackTop > 0) - { - stackTop--; - int index = stackIndex[stackTop]; - int exponent = stackExponent[stackTop]; - ulong product = stackProduct[stackTop]; - - if (index >= factorCount) - { - if (product != 1UL && product != order && candidateCount < limit) - { - ulong candidate = order / product; - if (candidate > 1UL && candidate < order) - { - candidates[candidateCount] = candidate; - candidateCount++; - } - } - - continue; - } - - int maxExponent = exponents[index]; - if (exponent > maxExponent) - { - continue; - } - - if (stackTop >= stackCapacity) - { - return candidateCount; - } - - stackIndex[stackTop] = index + 1; - stackExponent[stackTop] = 0; - stackProduct[stackTop] = product; - stackTop++; - - if (exponent == maxExponent) - { - continue; - } - - ulong primeFactor = factors[index]; - if (primeFactor == 0UL || product > order / primeFactor) - { - continue; - } - - if (stackTop >= stackCapacity) - { - return candidateCount; - } - - stackIndex[stackTop] = index; - stackExponent[stackTop] = exponent + 1; - stackProduct[stackTop] = product * primeFactor; - stackTop++; - } - - return candidateCount; - } - - private static void SortCandidatesKernel( - ulong prime, - ulong previousOrder, - bool hasPrevious, - ArrayView1D candidates, - int count) - { - for (int i = 1; i < count; i++) - { - ulong value = candidates[i]; - CandidateKey key = BuildCandidateKey(value, prime, previousOrder, hasPrevious); - int j = i - 1; - - while (j >= 0) - { - CandidateKey other = BuildCandidateKey(candidates[j], prime, previousOrder, hasPrevious); - if (other.CompareTo(key) <= 0) - { - break; - } - - candidates[j + 1] = candidates[j]; - j--; - } - - candidates[j + 1] = value; - } - } - - private static CandidateKey BuildCandidateKey(ulong value, ulong prime, ulong previousOrder, bool hasPrevious) - { - int group = GetGroup(value, prime); - if (group == 0) - { - return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); - } - - ulong reference = hasPrevious ? previousOrder : 0UL; - bool isGe = !hasPrevious || value >= reference; - int previousGroup = hasPrevious ? GetGroup(reference, prime) : 1; - int primary = ComputePrimary(group, isGe, previousGroup); - long secondary; - long tertiary; - - if (group == 3) - { - secondary = -(long)value; - tertiary = -(long)value; - } - else - { - ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; - secondary = (long)distance; - tertiary = (long)value; - } - - return new CandidateKey(primary, secondary, tertiary); - } - - private static int GetGroup(ulong value, ulong prime) - { - ulong threshold1 = prime >> 3; - if (value <= threshold1) - { - return 1; - } - - ulong threshold2 = prime >> 2; - if (value <= threshold2) - { - return 2; - } - - ulong threshold3 = (prime * 3UL) >> 3; - if (value <= threshold3) - { - return 3; - } - - return 0; - } - - private static int ComputePrimary(int group, bool isGe, int previousGroup) - { - int groupOffset; - switch (group) - { - case 1: - groupOffset = 0; - break; - case 2: - groupOffset = 2; - break; - case 3: - groupOffset = 4; - break; - default: - groupOffset = 6; - break; - } - - if (group == previousGroup) - { - if (group == 3) - { - return groupOffset + (isGe ? 0 : 3); - } - - return groupOffset + (isGe ? 0 : 1); - } - - return groupOffset + (isGe ? 0 : 1); - } - - private static bool TryConfirmCandidateKernel( - ulong prime, - ulong candidate, - in MontgomeryDivisorData divisor, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ArrayView1D compositeStack, - ArrayView1D statusOut, - ref int powUsed, - int powBudget) - { - int factorCount = FactorWithSmallPrimes(candidate, limit, primes, squares, primeCount, factors, exponents, out ulong remaining); - if (remaining != 1UL) - { - if (!TryFactorWithPollardKernel( - remaining, - limit, - primes, - squares, - primeCount, - factors, - exponents, - ref factorCount, - compositeStack, - statusOut)) - { - return false; - } - } - - SortFactors(factors, exponents, factorCount); - - for (int i = 0; i < factorCount; i++) - { - ulong primeFactor = factors[i]; - int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } - - ulong reduced = candidate; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (powBudget > 0 && powUsed >= powBudget) - { - return false; - } - - powUsed++; - if (Pow2EqualsOneKernel(reduced, divisor)) - { - return false; - } - } - } - - return true; - } - - private static bool TryFactorWithPollardKernel( - ulong initial, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ref int factorCount, - ArrayView1D compositeStack, - ArrayView1D statusOut) - { - if (initial <= 1UL) - { - return true; - } - - int stackCapacity = (int)compositeStack.Length; - if (stackCapacity <= 0) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - int stackTop = 0; - compositeStack[stackTop] = initial; - stackTop++; - - while (stackTop > 0) - { - stackTop--; - ulong composite = compositeStack[stackTop]; - if (composite <= 1UL) - { - continue; - } - - if (!PeelSmallPrimesKernel( - composite, - limit, - primes, - squares, - primeCount, - factors, - exponents, - ref factorCount, - out ulong remaining)) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - if (remaining == 1UL) - { - continue; - } - - ulong factor = PollardRhoKernel(remaining); - if (factor <= 1UL || factor == remaining) - { - if (!TryAppendFactorKernel(factors, exponents, ref factorCount, remaining, 1)) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - continue; - } - - ulong quotient = remaining / factor; - if (stackTop + 2 > stackCapacity) - { - statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow; - return false; - } - - compositeStack[stackTop] = factor; - compositeStack[stackTop + 1] = quotient; - stackTop += 2; - } - - return true; - } - - private static bool PeelSmallPrimesKernel( - ulong value, - uint limit, - ArrayView1D primes, - ArrayView1D squares, - int primeCount, - ArrayView1D factors, - ArrayView1D exponents, - ref int factorCount, - out ulong remaining) - { - ulong remainingLocal = value; - - for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate == 0U || primeCandidate > limit) - { - break; - } - - ulong square = squares[i]; - if (square != 0UL && square > remainingLocal) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remainingLocal % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remainingLocal /= primeValue; - exponent++; - } - while ((remainingLocal % primeValue) == 0UL); - - if (!TryAppendFactorKernel(factors, exponents, ref factorCount, primeValue, exponent)) - { - remaining = value; - return false; - } - } - - remaining = remainingLocal; - return true; - } - - private static bool TryAppendFactorKernel( - ArrayView1D factors, - ArrayView1D exponents, - ref int count, - ulong prime, - int exponent) - { - if (prime <= 1UL || exponent <= 0) - { - return true; - } - - for (int i = 0; i < count; i++) - { - if (factors[i] == prime) - { - exponents[i] += exponent; - return true; - } - } - - int capacity = (int)factors.Length; - if (count >= capacity) - { - return false; - } - - factors[count] = prime; - exponents[count] = exponent; - count++; - return true; - } - - private static ulong MulModKernel(ulong left, ulong right, ulong modulus) - { - GpuUInt128 product = new GpuUInt128(left); - return product.MulMod(right, modulus); - } - - private static ulong PollardRhoKernel(ulong value) - { - if ((value & 1UL) == 0UL) - { - return 2UL; - } - - ulong c = 1UL; - while (true) - { - ulong x = 2UL; - ulong y = 2UL; - ulong d = 1UL; - - while (d == 1UL) - { - x = AdvancePolynomialKernel(x, c, value); - y = AdvancePolynomialKernel(y, c, value); - y = AdvancePolynomialKernel(y, c, value); - - ulong diff = x > y ? x - y : y - x; - d = BinaryGcdKernel(diff, value); - } - - if (d == value) - { - c++; - if (c == 0UL) - { - c = 1UL; - } - - continue; - } - - return d; - } - } - - private static ulong AdvancePolynomialKernel(ulong x, ulong c, ulong modulus) - { - ulong squared = MulModKernel(x, x, modulus); - GpuUInt128 accumulator = new GpuUInt128(squared); - accumulator.AddMod(c, modulus); - return accumulator.Low; - } - - private static ulong BinaryGcdKernel(ulong a, ulong b) - { - if (a == 0UL) - { - return b; - } - - if (b == 0UL) - { - return a; - } - - int shift = BitOperations.TrailingZeroCount(a | b); - ulong aLocal = a >> BitOperations.TrailingZeroCount(a); - ulong bLocal = b; - - while (true) - { - bLocal >>= BitOperations.TrailingZeroCount(bLocal); - if (aLocal > bLocal) - { - ulong temp = aLocal; - aLocal = bLocal; - bLocal = temp; - } - - bLocal -= aLocal; - if (bLocal == 0UL) - { - return aLocal << shift; - } - } - } - - private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) - { - return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL; - } - - private static ulong CalculateByDoublingKernel(ulong prime) - { - ulong order = 1UL; - ulong value = 2UL % prime; - - while (value != 1UL) - { - value <<= 1; - if (value >= prime) - { - value -= prime; - } - - order++; - } - - return order; - } - private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan exponents, UInt128 prime, Span remainders) { if (exponents.Length == 0) @@ -1800,213 +658,26 @@ internal static partial class PrimeOrderGpuHeuristics { return Pow2ModKernelCache.GetOrAdd(accelerator, static accel => { - var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2ModKernel); + var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(PrimeOrderKernels.Pow2ModKernel); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); }); } - private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) - { - ulong exponent = exponents[index]; - remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false); - } - private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) { return Pow2ModKernelWideCache.GetOrAdd(accelerator, static accel => { - var loaded = accel.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(Pow2ModKernelWide); + var loaded = accel.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(PrimeOrderKernels.Pow2ModKernelWide); var kernel = KernelUtil.GetKernel(loaded); return kernel.CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); }); } - private static void Pow2ModKernelWide(Index1D index, ArrayView1D exponents, GpuUInt128 modulus, ArrayView1D remainders) - { - GpuUInt128 exponent = exponents[index]; - remainders[index] = Pow2ModKernelCore(exponent, modulus); - } - - private static GpuUInt128[] InitializeOddPowersTable(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) - { - GpuUInt128[] result = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount]; - result[0] = baseValue; - if (oddPowerCount == 1) - { - return result; - } - - // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. - baseValue.MulMod(baseValue, modulus); - - // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. - GpuUInt128 current = baseValue; - - // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. - for (int i = 1; i < oddPowerCount; i++) - { - current.MulMod(baseValue, modulus); - result[i] = current; - } - - return result; - } private static GpuUInt128 Pow2ModKernelCore(GpuUInt128 exponent, GpuUInt128 modulus) { - // This should never happen in production code. - // if (modulus == GpuUInt128.One) - // { - // return GpuUInt128.Zero; - // } - - // This should never happen in production code. - // if (exponent.IsZero) - // { - // return GpuUInt128.One; - // } - - GpuUInt128 baseValue = UInt128Numbers.TwoGpu; - - // This should never happen in production code - 2 should never be greater or equal modulus. - // if (baseValue.CompareTo(modulus) >= 0) - // { - // baseValue.Sub(modulus); - // } - - if (ShouldUseSingleBit(exponent)) - { - return Pow2MontgomeryModSingleBit(exponent, modulus, baseValue); - } - - int index = exponent.GetBitLength(); - int windowSize = GetWindowSize(index); - index--; - - int oddPowerCount = 1 << (windowSize - 1); - GpuUInt128[] oddPowers = InitializeOddPowersTable(baseValue, modulus, oddPowerCount); - GpuUInt128 result = GpuUInt128.One; - - int windowStart; - ulong windowValue; - while (index >= 0) - { - if (!IsBitSet(exponent, index)) - { - result.MulMod(result, modulus); - index--; - continue; - } - - windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } - - while (!IsBitSet(exponent, windowStart)) - { - windowStart++; - } - - // We're reusing oddPowerCount as windowBitCount here to lower registry pressure & avoid additional allocation - oddPowerCount = index - windowStart + 1; - index = windowStart - 1; - windowValue = ExtractWindowValue(exponent, windowStart, oddPowerCount); - - // We're reusing windowStart as square here to lower registry pressure & avoid additional allocation - for (windowStart = 0; windowStart < oddPowerCount; windowStart++) - { - result.MulMod(result, modulus); - } - - // We're reusing windowStart as tableIndex here to lower registry pressure & avoid additional allocation - windowStart = (int)((windowValue - 1UL) >> 1); - result.MulMod(oddPowers[windowStart], modulus); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool ShouldUseSingleBit(GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; - - private static int GetWindowSize(int bitLength) - { - if (bitLength <= Pow2WindowSizeBits) - { - return Math.Max(bitLength, 1); - } - - if (bitLength <= 23) - { - return 4; - } - - if (bitLength <= 79) - { - return 5; - } - - if (bitLength <= 239) - { - return 6; - } - - if (bitLength <= 671) - { - return 7; - } - - return Pow2WindowSizeBits; - } - - private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) - { - GpuUInt128 result = GpuUInt128.One; - - while (!exponent.IsZero) - { - if ((exponent.Low & 1UL) != 0UL) - { - result.MulMod(baseValue, modulus); - } - - exponent.ShiftRight(1); - if (exponent.IsZero) - { - break; - } - - // Reusing baseValue to store the squared base for the next iteration. - baseValue.MulMod(baseValue, modulus); - } - - return result; - } - - private static bool IsBitSet(GpuUInt128 value, int bitIndex) - { - if (bitIndex >= 64) - { - return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; - } - - return ((value.Low >> bitIndex) & 1UL) != 0UL; - } - - private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) - { - if (windowStart != 0) - { - exponent.ShiftRight(windowStart); - ulong mask = (1UL << windowBitCount) - 1UL; - return exponent.Low & mask; - } - - ulong directMask = (1UL << windowBitCount) - 1UL; - return exponent.Low & directMask; + return PrimeOrderKernels.Pow2ModKernelCore(exponent, modulus); } internal readonly record struct PrimeOrderGpuCapability(int ModulusBits, int ExponentBits) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 771fe28..3d364dc 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -847,7 +847,7 @@ public class MersenneDivisorCycles var kernel = accelerator.LoadAutoGroupedStreamKernel< Index1D, ArrayView1D, - ArrayView1D>(GpuDivisorCycleKernel); + ArrayView1D>(DivisorCycleKernels.GpuDivisorCycleKernel); using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true); @@ -941,93 +941,8 @@ public class MersenneDivisorCycles } - const byte ByteZero = 0; - const byte ByteOne = 1; - - // GPU kernel for divisor cycle calculation - static void GpuDivisorCycleKernel( - Index1D index, - ArrayView1D divisors, - ArrayView1D outCycles) - { - int i = index.X; - outCycles[i] = CalculateCycleLengthGpu(divisors[i]); - } - - // GPU-friendly version of cycle length calculation - /// - /// GPU-friendly cycle calculator that unrolls sixteen doubling steps; it wins the 8,388,607 benchmark and stays within ~2% of the octo loop at divisor 131,071. - /// - public static ulong CalculateCycleLengthGpu(ulong divisor) - { - if ((divisor & (divisor - 1UL)) == 0UL) - return 1UL; - - ulong order = 1UL; - ulong pow = 2UL; - - while (true) - { - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool GpuStep(ref ulong pow, ulong divisor, ref ulong order) - { - pow += pow; - if (pow >= divisor) - pow -= divisor; - - order++; - return pow == 1UL; - } + + public static ulong CalculateCycleLengthGpu(ulong divisor) => DivisorCycleKernels.CalculateCycleLengthGpu(divisor); public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false) { diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index ad53846..0a6bcc8 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -169,7 +169,7 @@ public sealed class PrimeTester(bool useInternal = false) { _accel = accelerator; // Compile once per accelerator and upload primes once. - Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(SmallPrimeSieveKernel); + Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); var primes = PrimesGenerator.SmallPrimes; DevicePrimes = accelerator.Allocate1D(primes.Length); @@ -266,56 +266,6 @@ public sealed class PrimeTester(bool useInternal = false) GpuKernelState.Clear(accelerator); } - // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise. - private static void SmallPrimeSieveKernel(Index1D index, ArrayView numbers, ArrayView smallPrimes, ArrayView results) - { - ulong n = numbers[index]; - if (n <= 3UL) - { - results[index] = (byte)(n >= 2UL ? 1 : 0); - return; - } - - if ((n & 1UL) == 0UL || (n > 5UL && (n % 5UL) == 0UL)) - { - // TODO: Replace the `% 5` branch with the GPU Mod5 helper once the sieve kernel can - // reuse the benchmarked bitmask reduction instead of modulo on every candidate. - results[index] = 0; - return; - } - - if (n.Mod10() == 1UL) - { - // Early reject special GCD heuristic with floor(log2 n) - ulong m = 63UL - (ulong)ILGPU.Algorithms.XMath.LeadingZeroCount(n); - if (BinaryGcdGpu(n, m) != 1UL) - { - results[index] = 0; - return; - } - } - - long len = smallPrimes.Length; - for (int i = 0; i < len; i++) - { - ulong p = smallPrimes[i]; - if (p * p > n) - { - break; - } - - if ((n % p) == 0UL) - { - // TODO: Route this modulo through the shared divisor-cycle cache once exposed to GPU kernels - // so batched sieves avoid per-prime `%` operations that profiling showed expensive. - results[index] = 0; - return; - } - } - - results[index] = 1; - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static bool SharesFactorWithMaxExponent(ulong n) { @@ -333,7 +283,7 @@ public sealed class PrimeTester(bool useInternal = false) // and divisor-cycle staging instead of allocating new device buffers per call. var gpu = GpuContextPool.Rent(); var accelerator = gpu.Accelerator; - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(SharesFactorKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel); int length = values.Length; var inputBuffer = accelerator.Allocate1D(length); @@ -384,45 +334,5 @@ public sealed class PrimeTester(bool useInternal = false) return u << shift; } - - private static void SharesFactorKernel(Index1D index, ArrayView numbers, ArrayView results) - { - ulong n = numbers[index]; - ulong m = 63UL - (ulong)XMath.LeadingZeroCount(n); - ulong gcd = BinaryGcdGpu(n, m); - results[index] = gcd == 1UL ? (byte)0 : (byte)1; - } - - private static ulong BinaryGcdGpu(ulong u, ulong v) - { - // TODO: Replace this inline GPU binary GCD with the kernel extracted from - // GpuUInt128BinaryGcdBenchmarks via GpuKernelPool so device callers reuse the - // fully unrolled ladder instead of this branchy fallback. - if (u == 0UL) - { - return v; - } - - if (v == 0UL) - { - return u; - } - - int shift = XMath.TrailingZeroCount(u | v); - u >>= XMath.TrailingZeroCount(u); - - do - { - v >>= XMath.TrailingZeroCount(v); - if (u > v) - { - (u, v) = (v, u); - } - - v -= u; - } - while (v != 0UL); - - return u << shift; } } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 9351b65..fe97565 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -536,7 +536,7 @@ public static class ULongExtensions var kernel = KernelCache.GetOrAdd(accelerator, static accel => { - var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernel); + var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernel); var launcher = KernelUtil.GetKernel(loaded); return launcher.CreateLauncherDelegate, MontgomeryDivisorData, byte, ArrayView1D>>(); }); @@ -573,11 +573,6 @@ public static class ULongExtensions } } - private static void Pow2MontgomeryKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, byte keepMontgomery, ArrayView1D results) - { - results[index] = Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0); - } - } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) diff --git a/tools/GpuNumericTypeProbe.csx b/tools/GpuNumericTypeProbe.csx index fe3d284..7e97a76 100644 --- a/tools/GpuNumericTypeProbe.csx +++ b/tools/GpuNumericTypeProbe.csx @@ -1,6 +1,7 @@ #r "nuget: ILGPU, 1.5.3" #r "nuget: PeterO.Numbers, 1.8.1" #r "nuget: Open.Numeric.Primes, 4.0.4" +#load "../PerfectNumbers.Core/Gpu/Kernels/NumericProbeKernels.cs" using System; using System.Numerics; @@ -9,6 +10,7 @@ using ILGPU.Runtime; using ILGPU.Runtime.CPU; using PeterO.Numbers; using Open.Numeric.Primes; +using PerfectNumbers.Core.Gpu; var context = Context.CreateDefault(); var accelerator = context.CreateCPUAccelerator(0); @@ -53,7 +55,7 @@ static ProbeResult ProbeBigInteger(Accelerator accelerator) { try { - var kernel = accelerator.LoadAutoGroupedStreamKernel>(BigIntegerKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel>(NumericProbeKernels.BigIntegerKernel); using var buffer = accelerator.Allocate1D(1); var extent = new Index1D((int)buffer.Length); kernel(extent, buffer.View); @@ -72,7 +74,7 @@ static ProbeResult ProbeEInteger(Accelerator accelerator) { try { - var kernel = accelerator.LoadAutoGroupedStreamKernel>(EIntegerKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel>(NumericProbeKernels.EIntegerKernel); using var buffer = accelerator.Allocate1D(1); kernel(new Index1D((int)buffer.Length), buffer.View); accelerator.Synchronize(); @@ -92,7 +94,7 @@ static ProbeResult ProbeERational(Accelerator accelerator) try { - var kernel = accelerator.LoadAutoGroupedStreamKernel>(ERationalKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel>(NumericProbeKernels.ERationalKernel); using var buffer = accelerator.Allocate1D(1); kernel(new Index1D((int)buffer.Length), buffer.View); accelerator.Synchronize(); @@ -119,7 +121,7 @@ static ProbeResult ProbeOpenNumericPrimesIsPrime(Accelerator accelerator) { try { - var kernel = accelerator.LoadAutoGroupedStreamKernel>(OpenNumericPrimesIsPrimeKernel); + var kernel = accelerator.LoadAutoGroupedStreamKernel>(NumericProbeKernels.OpenNumericPrimesIsPrimeKernel); using var buffer = accelerator.Allocate1D(4); kernel(new Index1D((int)buffer.Length), buffer.View); accelerator.Synchronize(); @@ -133,47 +135,6 @@ static ProbeResult ProbeOpenNumericPrimesIsPrime(Accelerator accelerator) } } -static void OpenNumericPrimesIsPrimeKernel(Index1D index, ArrayView1D data) -{ - ulong value = (ulong)(index + 2); - bool isPrime = Prime.Numbers.IsPrime(value); - data[index] = isPrime ? 1 : 0; -} - -static void BigIntegerKernel(Index1D index, ArrayView1D data) -{ - BigInteger value = new BigInteger(data[index]); - value += BigInteger.One; - value -= BigInteger.One; - value *= new BigInteger(2); - value /= new BigInteger(3); - value %= new BigInteger(5); - int comparison = value.CompareTo(BigInteger.Zero); - data[index] = comparison; -} - -static void EIntegerKernel(Index1D index, ArrayView1D data) -{ - EInteger value = EInteger.FromInt32(data[index]); - value = value + EInteger.One; - value = value - EInteger.One; - value = value * EInteger.FromInt32(2); - value = value.Divide(EInteger.FromInt32(3)); - value = value.Remainder(EInteger.FromInt32(5)); - int comparison = value.CompareTo(EInteger.Zero); - data[index] = comparison; -} - -static void ERationalKernel(Index1D index, ArrayView1D data) -{ - ERational value = ERational.FromInt32(data[index]); - ERational adjusted = (value + ERational.FromInt32(1)) - ERational.FromInt32(1); - ERational product = adjusted * ERational.FromInt32(2); - ERational quotient = product / ERational.FromInt32(3); - int comparison = quotient.CompareTo(ERational.FromInt32(0)); - data[index] = comparison; -} - static bool HasERationalRemainder() { foreach (var method in typeof(ERational).GetMethods())