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 1/2] 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()) From 473d4efcdfa8457017e7d865e0378f792ea61f47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sat, 18 Oct 2025 17:43:24 +0200 Subject: [PATCH 2/2] Clarify numeric helper placement guidance --- AGENTS.md | 290 ++++++++++++++++++++++++------------------------------ 1 file changed, 127 insertions(+), 163 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 35271f7..2caac30 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -1,191 +1,155 @@ # Repository Guidelines -This repository contains only one solution: +- Only solution: `EvenPerfectScanner.sln` (console app scanning numeric ranges). +- Keep changes confined to the relevant solution/project unless instructed otherwise. +- Assume each solution has or will have its own tests and build configuration; avoid cross-project edits unless required. -- `EvenPerfectScanner.sln` – a console application that scans numeric ranges. +## Build and Test -Changes should be isolated to the relevant solution and project. When working on one -project, avoid modifying the other unless explicitly required. Each solution -has or will have its own test suite and build configuration. +- Build: `dotnet build EvenPerfectScanner.sln`. +- Run only tests related to the modified solution. +- Run tests without `--no-build`. +- Write tests with xUnit + FluentAssertions. +- Confirm test runtimes are reasonable; skip suites expected to run excessively long. +- Always apply timeout guards and `Category=Fast` filters when running tests. -To build `EvenPerfectScanner` run: +## File Editing -```bash -dotnet build EvenPerfectScanner.sln -``` +- Do not use `apply_patch.cs` until TODOs lift the ban; edit files with shell tools and stage via git. +- Intended usage (once re-enabled): `dotnet script apply_patch.cs -- [target-directory]`. +- Expect the script to run in local/remote Codex environments on Windows and Linux. +- Supported options (when allowed): + - `--check` + - `--ignore-whitespace` + - `--ignore-eol` + - `--ignore-line-endings` + - `--target ` + - `--profile` + - `--replace-line ` (use `_` for expected to skip validation) +- Build patches in standard `git diff` format: + - Include `diff --git`, `---`, `+++`, and `@@` headers with numeric ranges. + - Exclude `***` or `*** Begin/End Patch` wrappers. +- Run `--check` before applying patches; apply only after a successful check. +- Keep patches small, precise, and regenerated after each edit; confirm target lines with `Get-Content` or `nl -ba`. +- Record exact start lines for each hunk and double-check with `nl -ba`. +- Prefer `--replace-line` for single-line edits after capturing the target line number. +- Regenerate diffs after successful patch applications and whenever files shift. +- Use `--ignore-whitespace` / `--ignore-line-endings` only when validation noise blocks `--check`. +- Maintain unique patch filenames and delete applied ones. +- Investigate failed hunks immediately and rebuild diffs instead of retrying stale context. +- Prefer single-purpose patches. +- Use `nl -ba` to show line numbers when needed. +- Use four spaces for code block indentation. +- Create patches from the repo root with relative paths. +- Do not apply changes manually or via other scripts unless modifying `apply_patch.cs`. +- If `apply_patch.cs` has issues or lacks features, fix or enhance it and continue using it. -Run only the tests related to the solution you modify. Always run tests without --no-build parameter to avoid run failures due to missing assembly files. Tests should be always written using XUnit with FluentAssertions. Before executing any tests, make sure their running time is reasonable; skip test suites that are expected to take excessively long. When running tests, always provide both an explicit timeout guard and the Category=Fast filter so that only the fast suite executes. +## Test Execution Time Policy -## Making Changes in the Repository -To make any changes to the files in the repository, you should always utilize a script created for it to overcome environment / scripting issues. The script is called `apply_patch.cs` and has the following syntax: +- Estimate runtimes before launching tests or samples. +- Do not run tests expected to exceed 2 minutes unless explicitly authorized. +- Never impose the two-minute guard on benchmarks or profiling runs; let them finish unless obviously stuck. +- Keep timeouts external to the test code. +- Use targeted filters such as `--filter FullyQualifiedName~` or `--filter "FullyQualifiedName=Namespace.Class.Method"`. +- Keep GPU-heavy tests minimal; skip kernels needing >2 minutes for compile/warmup. +- In local Codex tasks, run only fast, newly added tests and list recommended existing tests instead. -```bash -dotnet script apply_patch.cs -- [target-directory] -``` - -> **Temporary workflow notice:** Until the pending apply_patch.cs performance and correctness TODOs are complete, do not invoke the script in this repository. Edit files directly with shell tools (sed, tee, cat, python, etc.) and stage the results via git. Re-enable the helper only after the TODO tracker explicitly calls out that the ban has been lifted. - -It supports the standard `git diff` patch files as input. You should be able to run it both in the local and remote `Codex` environments. You should expect it to work both in `Windows` and `Linux` environments. Make sure that you include `diff --git` lines, headers and file prefixes in the patch files to avoid applying issues. - -`apply_patch.cs` supports several parameters, depending on your needs: -* `--check` - allows checking if the patch will work without making changes to the files -* `--ignore-whitespace` -* `--ignore-eol` -* `--ignore-line-endings` -* `--target` - allows specifying target root directory for the patch -* `--profile` - capture read/parse/apply timings for each patch file. -* `--replace-line ` - replace a single line by absolute number (use `_` as the expected text to skip validation) - -* Build patches in standard `git diff` format. - * It supports `---`, `+++`, and `@@`. - * It doesn't support `***` or `*** Begin/End Patch` wrappers. - * It requires `diff --git` at the beginning. - * Patch header needs the actual line numbers. - * When hand-writing a diff, always include the numeric ranges in every `@@ -start,count +start,count @@` header; missing or zero counts will cause `--check` to fail immediately. -* Always check if the patch will work by using `--check` parameter, without making changes to the files, which would force patch changes. Try applying the patch only when `--check` passed. -* Build patches in small, line-precise hunks; confirm target lines with `Get-Content` or `nl -ba` before writing the diff to avoid context mismatches. -* Capture the exact start lines for each hunk while drafting the `@@ -start,length +start,length @@` header; double-check them with `nl -ba` so the first verification succeeds. -* For one-off line edits, prefer the `--replace-line ` option. Pass `_` as the expected text to skip validation when necessary. -* Before calling `--replace-line`, capture the target number with `nl -ba` (or an equivalent command) so the replacement lands exactly once without retrying. -* After every successful patch application, regenerate a fresh diff from the updated files before preparing the next patch so the hunk headers and context reflect the latest line numbers. -* If the patch script flags repeated mismatches, regenerate a fresh diff for the updated file instead of tweaking stale context blocks. -* When `--check` reveals that only whitespace or EOL normalization differs, rerun it with `--ignore-whitespace` and/or `--ignore-line-endings` so validation concentrates on the meaningful edits. -* Regenerate the diff whenever the file shifts—never reuse an old hunk after other edits. -* Keep patch filenames unique and delete applied files to prevent accidental re-use. -* When a hunk fails, inspect the around-lines immediately and adjust rather than retrying the same diff. -* Prefer single-purpose patches (one logical change per file) so rollbacks or fixes stay focused. -* Use `nl -ba` to show line numbers, when needed. -* Always use 4 spaces for indentation of code blocks. - -If you want to make any changes, create a patch file and run the script to apply it. Always run it from the root directory and reference any modified files with relative paths respectively in the patch. - -Do **not** attempt to apply the changes manually, or with `Python`, or with `PowerShell` scripts, or `cat` or any other way, unless you apply changes to `apply_patch.cs` script itself. If you identify any issues with `apply_patch.cs` script, you should propose fixes to resolve them and continue using the script after updates. If you spot any missing, but required features, propose enhancements and continue using the script after updates, too. - -### Test execution time policy - -- Always estimate runtime before launching any tests or long-running samples. If a test or ad-hoc run is likely to exceed 2 minutes, do not run it unless the user explicitly requests a longer run and provides justification. -- Never apply the two-minute guard to benchmark or profiling runs; allow them to finish and capture the summary unless the console log shows the job is clearly stuck. -- Do not hardcode timeouts inside test code. Enforce time limits from the runner/shell layer so developers can run tests locally without artificial constraints. -- Prefer targeted filters over full suites to keep runs short: `--filter FullyQualifiedName~` or `--filter "FullyQualifiedName=Namespace.Class.Method"`. -- Keep GPU-heavy tests minimal and scoped. If a kernel compile or warmup is expected to take >2 minutes, skip running it during regular iterations. -- When working on a task in a local Codex environment, do not run any existing unit tests, except for fast, new, added tests, unless you're instructed otherwise. In the task summary list the names of the recommended unit tests to validate the changes, instead. - -#### How to run tests within the 2-minute limit +### 120-second Guard Examples - Windows (PowerShell): - - Run a single class or test with a 120s guard: - - `$proj = ".csproj"` - - `$filter = 'FullyQualifiedName~SomeTests' # or FullyQualifiedName=Namespace.Class.Method` - - `$job = Start-Job -ScriptBlock { param($p,$f) dotnet test $p -c Debug --filter $f } -ArgumentList $proj,$filter` - - `if (-not (Wait-Job $job -Timeout 120)) { Stop-Job $job; throw 'Test run exceeded 120s timeout.' }` - - `Receive-Job $job` - + - `$proj = ".csproj"` + - `$filter = 'FullyQualifiedName~SomeTests'` + - `$job = Start-Job -ScriptBlock { param($p,$f) dotnet test $p -c Debug --filter $f } -ArgumentList $proj,$filter` + - `if (-not (Wait-Job $job -Timeout 120)) { Stop-Job $job; throw 'Test run exceeded 120s timeout.' }` + - `Receive-Job $job` - Linux/macOS (bash): - - Use GNU timeout (install coreutils if needed on macOS): - - `timeout 120s dotnet test .csproj -c Debug --filter "Category=Fast&FullyQualifiedName~SomeTests"` - - Always specify `timeout 120s` (or a shorter limit when appropriate) and include `Category=Fast` in every filter. - - Note: exit code 124 indicates timeout; handle accordingly in scripts. + - `timeout 120s dotnet test .csproj -c Debug --filter "Category=Fast&FullyQualifiedName~SomeTests"` + - Always specify `timeout 120s` (or shorter) and include `Category=Fast` in the filter. + - Treat exit code 124 as a timeout. +- Codex/agent runs: wrap `dotnet test` commands in shell guards above. +- Filter examples: `--filter "FullyQualifiedName~PrimeTesterGcdTests"`, `--filter "FullyQualifiedName=EvenPerfectBitScanner.Tests.ResultsFileNameTests.BuildResultsFileName_contains_the_range"`. -- Codex/agent runs: when invoked as separate steps where inline time control isn’t available, wrap calls through the shell using the guards above rather than calling `dotnet test` directly from the tool. +## Coding Standards -Examples: -- Quick GCD check: `--filter "FullyQualifiedName~PrimeTesterGcdTests"` -- Single test method: `--filter "FullyQualifiedName=EvenPerfectBitScanner.Tests.ResultsFileNameTests.BuildResultsFileName_contains_reduction_and_devices"` +- Use four-space indentation; increase one level per nested block. +- Indent wrapped parameters/expressions by two indentation units. +- Reuse existing variables when prior values are obsolete; call out meaning changes in comments to keep intent clear and register pressure low. +- Avoid try/finally blocks solely for resource release unless omitting finally harms clarity or maintainability; ensure every path releases resources. +- Limit the divisor cycle cache to one in-memory block plus the on-disk snapshot; recompute on misses, discard transient results, and avoid background generation or extra caches. +- Write code, comments, commit messages, branch names, and PR descriptions in American English. +- Preserve relevant TODO notes and explanatory comments unless invalid. +- Prefer structs/classes with public fields over auto-properties for hot-path data holders. +- Pass mutable structs as parameters with the `in` modifier whenever the call site and usage allow it (skip when constraints such as lambda captures require a copy). +- Omit the `in` modifier when passing immutable structs as parameters. +- Remove or extend trivial wrapper methods that only call another method. -Always implement the functionality to achieve the highest possible performance, even if that means allowing for some code duplication. E.g. use List instead of IList, IEnumerable etc.; calculate modulo "%" result in the most efficient way, use aggresive inlining when appropriate and so on. Always take arrays from pools, when it's desired and possible depending on the required array size. +## Performance Checklist -Don't simplify the functionality and / or algorithms to behave differently, unless you're resolving a bug that you've found. Always attempt preservint the functionality, is it's designed in iterations and there is a reason for specific behaviors / implementations. +- Favor dedicated kernels over flag-controlled/parameterized ones; duplicate code when it improves speed. +- Remove `if` statements inside hot loops when possible. +- Replace runtime flags in hot paths with dedicated method variants. +- Reuse variables instead of allocating new ones; document meaning changes. +- Extract functionality into separate classes and place reusable parts in `*.Core` projects. +- Rent arrays from `ArrayPool` when stack allocation is impractical and return them after use. +- Provide cancellation-aware and cancellation-free variants of hot methods; handle cancellation in outer scopes. +- Prefer `for` loops over `foreach` to avoid enumerator allocations. +- Declare loop control variables outside loops and reuse them. +- Cache arrays, static/read-only arrays, and collection counts in locals before loops. +- Declare variables together outside loops when feasible and reuse them. +- Avoid repeated casting within loops by precomputing or using wider variables. -The only .Net SDK available in Codex environment is .Net 8.0. It'll be periodically updated to newer versions as they become available. Specify target framework to the one currently available, but always target all .csproj files against all .Net versions currently supported by MS, using the latest official release. .Net 6.0 is no longer supported, so projects should target only .Net 8.0 and .Net 9.0, with 8+ being the default. +## GPU Notes -Always enter blank line after closing "}" parenthesis, unless there is a consecutive closing "}" bracket. IF's, loops, in practice all code blocks should be put in their "{ (...) }" blocks. Even if there is only a single line in such a code block. Always use the latest officially released language features. All .csproj files should define nullable = enabled and we should always use "?" and/or "!" where appropriate. +- ILGPU 1.5.3 checks: `%` and integer division supported for 32-bit/64-bit integers; `UInt128` remains unsupported inside kernels. Use the `GpuUInt128` struct for GPU-side emulation and extend it whenever required arithmetic helpers are missing. +- Arithmetic in GPU kernels covers only 32-bit and 64-bit primitives; rely on `GpuUInt128` together with its read-only counterpart `ReadOnlyGpuUInt128` to emulate 128-bit behavior. +- CPU code should default to `UInt128` instead of `GpuUInt128` unless the same calculation executes millions or billions of times, in which case prefer mutable `GpuUInt128` when benchmarks show no material regression; benchmark summaries and recommendations live in `Benchmarks.md`. +- GPU kernels should favor mutable `GpuUInt128` and `ReadOnlyGpuUInt128` when the necessary operations exist on the accelerator and benchmarks do not show a dramatic advantage for `UInt128`. +- Perform conversions from `GpuUInt128` or other structs to `ReadOnlyGpuUInt128` outside hot loops whenever possible; push the conversion toward the top of the call chain to avoid repeated costs. +- For each numeric type, review available helpers in `UIntExtensions`, `ULongExtensions`, `UInt128Extensions`, and `GpuUInt128`; add new operations to the extension class or struct that matches the type. +- Arrays/vectors of `bool` unsupported; use `byte`. +- `UInt128` multiplication, division, bit shifts unsupported in kernels; implement manually. +- GPU context pooling supported and recommended. +- Numeric probe results (ILGPU 1.5.3, CPU accelerator): + - `System.Numerics.BigInteger`: unsupported (`InternalCompilerException`). + - `PeterO.Numbers.EInteger`: unsupported (`InternalCompilerException`). + - `PeterO.Numbers.ERational`: unsupported (`InternalCompilerException`, lacks remainder API). -Never ever run tests marked with trait Category = Slow. Always invoke `dotnet test` (or equivalent wrappers) with both an explicit timeout and the `Category=Fast` filter; never run tests without specifying these arguments. If there are no categories assigned yet, categorize them before running or skip the execution. +## Patch Preparation Reminders -There is only one agents file - this one. Don't search for another one. +- Manually review every patch before `--check`; remove unintended whitespace noise. +- Capture context with `nl -ba` before composing each hunk; include ≥3 context lines when available. +- Normalize trailing whitespace early or rely on ignore flags temporarily, then rerun without ignores to confirm the final diff. +- Regenerate patches from scratch after working tree changes; avoid editing old diffs. +- Note that the parser ignores blank metadata lines between headers/hunks; keep headers compact unless tools add blanks. -Always use indentation with four space characters in source code files, even if a file currently uses tabs or appears to rely on eight-space indentation. Increase indentation level by exactly one indentation unit (four spaces) per nested block. When parameters or long expressions wrap, indent the continuation with two indentation units for clarity. -Explicitly prefer reusing existing variables once their previous values are no longer required instead of declaring new locals of the same type. Call out the variable reuse in comments when its meaning changes so the intent stays clear and register pressure remains minimal. -Do not introduce try/finally blocks solely to release pooled resources. Restructure the code so that every return path releases the resources without relying on finally unless the control flow would become unclear or unmaintainable. +## apply_patch.cs Profiling Snapshot (2025-10-04) -Keep the divisor cycle cache limited to a single block in memory (only the on-disk snapshot). When a lookup misses, compute the required cycle on the configured device and discard the transient result instead of scheduling background generation or caching additional blocks. - -All code, comments, commit messages, branch names and PR descriptions must be written in American English. - -- When modifying existing files, preserve all relevant TODO notes and explanatory comments unless they are no longer valid. -- Prefer defining data holders as structs or classes with public fields over auto-properties when performance on the hot path is a concern. -- Do not add or keep methods that simply wrap another method call without adding any logic; inline those calls or extend the wrapper with meaningful work. - -## Performance implementation checklist - -- Prefer dedicated kernels over flag-controlled or parameterized kernels. Avoid conditional branches in kernels and tight loops by creating separate kernels when necessary. -- Remove `if` statements inside hot loops whenever possible. Duplicated code in additional kernels or methods is acceptable when it improves speed. -- Replace runtime flags in hot paths with dedicated method variants to eliminate conditional branches. -- Reuse existing variables of the same type instead of allocating new ones if some become unused. Add comments whenever a variable's meaning changes. -- Extract functionality into separate classes instead of growing existing classes. Place reusable components into `*.Core` projects. -- Rent arrays from `ArrayPool` when stack allocation is impractical or less efficient, and return them after use. -- Provide cancellation-aware and cancellation-free variants of hot methods so cancellation tokens are not checked in every loop iteration; handle cancellation in outer scopes. -- Prefer `for` loops over `foreach` to reduce enumerator allocations. -- Declare loop control variables outside the loop and reuse them across loops when possible. -- Cache static or read-only arrays in local variables before looping to avoid repeated class field lookups. -- Store collection `Count` or array `Length` in a local variable before the loop and use the cached value. -- Declare variables together outside loops when feasible and reuse them instead of creating new ones. -- Avoid repeated casting within loops by precomputing or using wider variables when necessary. - -## GPU-specific notes - -The following checks were executed with ILGPU 1.5.3 on CPU and OpenCL accelerators: - -- [x] `%` operator (kernels): supported for 32-bit and 64-bit integer types; not supported for `UInt128` inside kernels. -- [x] Integer division (kernels): supported for 32-bit and 64-bit types; not supported for `UInt128` inside kernels. -- [x] Arrays or vectors of `bool`: unsupported (`NotSupportedException`); use `byte` instead. -- [x] `UInt128` multiplication: unsupported in kernels (`AggregateException/NullReferenceException`); manual implementation required. -- [x] `UInt128` division: unsupported in kernels; manual implementation required. -- [x] `UInt128` bit shifts: attempts failed in kernels; avoid inside device code. -- [x] GPU context pooling: supported and recommended (validated in a separate program). - -Additional numeric-type probe results (`tools/GpuNumericTypeProbe.csx`, ILGPU 1.5.3 on the CPU accelerator): - -- `System.Numerics.BigInteger`: unsupported. Kernels that exercise addition, subtraction, multiplication, division, modulo, and `CompareTo` all fail with `InternalCompilerException`, indicating fundamental codegen gaps for the type. -- `PeterO.Numbers.EInteger`: unsupported. Attempts to compile kernels that call addition, subtraction, multiplication, division, remainder, and `CompareTo` throw `InternalCompilerException`, so none of the arithmetic helpers are currently usable on the device. -- `PeterO.Numbers.ERational`: unsupported. Kernels that call addition, subtraction, multiplication, division, and `CompareTo` trigger `InternalCompilerException`; the type also lacks a remainder/modulo API for device usage. - -Patch-preparation checklist for reliable `--check` runs: - -- Manually review every patch file before running `--check`; confirm the hunks reflect only the intended edits and eliminate accidental whitespace-only noise ahead of time. -- Always capture context with `nl -ba` before composing each hunk so the `@@` header line numbers match the target file on the first try. -- Ensure every hunk carries at least three unchanged context lines above and below the edits whenever the surrounding file provides them; regenerate the diff if the file changed since the previous attempt so the context stays accurate. -- Trim or normalize trailing whitespace while drafting the patch, or rely on `--ignore-whitespace` during verification when formatting noise is unavoidable, then re-run without ignores to confirm the final diff matches exactly. -- Prefer regenerating the patch from scratch whenever the working tree changes instead of editing previous diffs; stale context is the most common source of mismatches. - -- The patch parser now ignores blank metadata lines between `diff --git`, `---`, and `+++` headers and skips stray blank separators between hunks, so you no longer need to insert dummy empty lines to satisfy the checker. Keep the header compact unless your diff tool emits those blank lines naturally. - -### apply_patch.cs profiling snapshot (2025-10-04) - -- Latest timing check: `dotnet script apply_patch.cs -- --profile --check /tmp/agents_todo_update.diff`. -- Patch characteristics: 2 files, 4 hunks, 31 lines, 3468 bytes. +- `dotnet script apply_patch.cs -- --profile --check /tmp/agents_todo_update.diff`. +- Sample patch: 2 files, 4 hunks, 31 lines, 3468 bytes. - Timings: read 0.45 ms, parse 2.96 ms, apply 2.96 ms. -- Focus next on streaming the patch reader and pooled buffer reuse to reduce the parse APPLY bottleneck highlighted above. +- Next focus: stream the patch reader and reuse pooled buffers to reduce parse/apply cost. -Mersenne divisor size estimates: +## Mersenne Divisor Guidance -- Any prime divisor `q` of a Mersenne number `2^p - 1` must satisfy `q ≡ 1 (mod 2p)`. The minimal `q` therefore appears at `q = 2p + 1` when the multiplier `k = 1` yields a prime. -- For candidate exponents `p` in the `[138,000,000, 200,000,000]` range, the minimal compliant divisor lands below `2 · 200,000,000 + 1 = 400,000,001`, which fits in 29 bits (`2^29 = 536,870,912`). -- Extending the search up to `p = 600,000,000` pushes the minimal divisor to `1,200,000,001`, requiring 31 bits (`2^31 = 2,147,483,648`). -- Larger `k` values (and thus larger `q`) are possible—up to nearly `2^p - 1`—so the heuristics and storage should plan for wider integers when deep trial factoring is necessary, but the noted bit widths cover the first admissible divisors for these exponent bands. +- Prime divisor `q` of `2^p - 1` satisfies `q ≡ 1 (mod 2p)` with minimal `q = 2p + 1` for `k = 1`. +- Target for at least `p` in `[138,000,000, 200,000,000]`, minimal `q < 400,000,001` (29 bits). +- Target for at least `p = 600,000,000`, minimal `q ≤ 1,200,000,001` (31 bits). +- Larger `k` values can grow `q` up to nearly `2^p - 1`; plan heuristics/storage for wider integers even though the noted bit widths cover the first admissible divisors. +- Worst-case divisor width: `p - 1` bits (`⌊(2^p - 1) / 2⌋ = 2^{p-1} - 1`). +- `k_max ≈ 2^{p-2} / p` for large `p`. -- Because every proper divisor of `2^p - 1` is at most `(2^p - 1) / 2`, the worst-case bit-width for any admissible `q` is `p - 1` bits (since `⌊(2^p - 1) / 2⌋ = 2^{p-1} - 1`). Consequently we need roughly `199,999,999` bits when `p ≤ 200,000,000` and `599,999,999` bits when `p ≤ 600,000,000`. In terms of the `q = 2kp + 1` form, this translates to `k_max = ⌊((2^p - 1) / 2 - 1) / (2p)⌋`, which is approximately `2^{p-2} / p` for large `p`. +## Practical Kernel Guidance -Practical guidance: -- Use native `%` and `/` in kernels when working with 32-bit or 64-bit integers if it is the most efficient option for the target backend. -- For `UInt128` in kernels, avoid `%`/`/`; implement reductions via Montgomery/Barrett, or use Mersenne-specific folding for mod `2^p-1`. -- Prefer Mersenne folding over general modulo for `2^p-1` moduli to remove branches in hot loops. +- Use native `%` and `/` for 32-bit/64-bit integers when efficient. +- For `UInt128` in kernels, avoid `%`/`/`; prefer Montgomery/Barrett or Mersenne folding. +- Prefer Mersenne folding for `2^p - 1` moduli to remove branches in hot loops. -- Heuristic multiply-shift benchmarks (2025-10-04): - - Run `dotnet run -c Release --project tools/HeuristicMultiplyShiftTiming/HeuristicMultiplyShiftTiming.csproj` to compare the `ULongExtensions` helpers without the BenchmarkDotNet harness. - - Timings below report the average nanoseconds per call over 200,000,000 iterations; all variants produced identical checksums, confirming consistent arithmetic. - - `NearOverflowShift`: `MultiplyShiftRight` ≈ 2.37 ns, shift-first ≈ 4.71 ns, naive `(value * multiplier) >> shift` ≈ 0.85 ns. - - `HalfRange`: `MultiplyShiftRight` ≈ 2.11 ns, shift-first ≈ 4.15 ns, naive ≈ 0.70 ns. - - `MixedBits`: `MultiplyShiftRight` ≈ 1.43 ns, shift-first ≈ 2.83 ns, naive ≈ 0.60 ns. - - `Sparse`: `MultiplyShiftRight` ≈ 1.57 ns, shift-first ≈ 2.84 ns, naive ≈ 0.59 ns. - - The shift-first rewrite roughly doubles the per-call cost on CPU, but it stays within 64-bit arithmetic and is the preferred fallback for GPU heuristics when avoiding `UInt128` overflow is mandatory. +## Heuristic Multiply-Shift Benchmarks (2025-10-04) + +- Run `dotnet run -c Release --project tools/HeuristicMultiplyShiftTiming/HeuristicMultiplyShiftTiming.csproj` to compare `ULongExtensions` helpers without the BenchmarkDotNet harness. +- Average ns per call over 200,000,000 iterations (identical checksums): + - `NearOverflowShift`: `MultiplyShiftRight` ≈ 2.37 ns, shift-first ≈ 4.71 ns, naive ≈ 0.85 ns. + - `HalfRange`: `MultiplyShiftRight` ≈ 2.11 ns, shift-first ≈ 4.15 ns, naive ≈ 0.70 ns. + - `MixedBits`: `MultiplyShiftRight` ≈ 1.43 ns, shift-first ≈ 2.83 ns, naive ≈ 0.60 ns. + - `Sparse`: `MultiplyShiftRight` ≈ 1.57 ns, shift-first ≈ 2.84 ns, naive ≈ 0.59 ns. +- Shift-first doubles CPU cost but stays within 64-bit arithmetic and is the GPU fallback when avoiding `UInt128` overflow.