diff --git a/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs index e7c0c7a..73ce540 100644 --- a/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs @@ -1,5 +1,6 @@ using System.Runtime.CompilerServices; using BenchmarkDotNet.Attributes; +using PerfectNumbers.Core; using BenchmarkDotNet.Jobs; namespace EvenPerfectBitScanner.Benchmarks; @@ -183,10 +184,11 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_Current(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; @@ -208,10 +210,11 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_DoWhile(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; @@ -236,10 +239,11 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_DoubleSubtract(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; @@ -265,10 +269,11 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_UInt128Modulo(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; @@ -287,22 +292,23 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_UnrolledPair(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; while (true) { - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } @@ -312,32 +318,33 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_UnrolledQuad(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; while (true) { - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } @@ -347,52 +354,53 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_UnrolledOct(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; while (true) { - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } @@ -402,101 +410,101 @@ public class MersenneDivisorCycleLengthGpuBenchmarks [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong CalculateCycleLengthGpu_UnrolledHex(ulong divisor) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - return 1UL; - } + // EvenPerfectBitScanner never routes power-of-two divisors here; keep the guard disabled outside tests and benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } ulong order = 1UL; ulong pow = 2UL; while (true) { - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - if (Step(ref pow, divisor, ref order)) + if (AdvanceMersenneDivisorCycleStep(ref pow, divisor, ref order)) { return order; } - } } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool Step(ref ulong pow, ulong divisor, ref ulong order) + private static bool AdvanceMersenneDivisorCycleStep(ref ulong pow, ulong divisor, ref ulong order) { pow += pow; if (pow >= divisor) @@ -507,4 +515,5 @@ public class MersenneDivisorCycleLengthGpuBenchmarks order++; return pow == 1UL; } + } diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 3d364dc..ae214cf 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -206,9 +206,6 @@ public class MersenneDivisorCycles } } - // TODO: Replace this naive doubling fallback with the unrolled-hex generator from - // MersenneDivisorCycleLengthGpuBenchmarks so large divisors avoid the millions of - // iterations measured in the scalar loop. // TODO: Route this miss to an ephemeral single-cycle computation on the device selected by the current // configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot // so we respect the large-p memory limits and avoid extra cache locks. @@ -987,25 +984,11 @@ public class MersenneDivisorCycles private static ulong CalculateCycleLengthFallback(ulong divisor) { - // Otherwise, find order of 2 mod divisor - // TODO: Switch this scalar fallback to the unrolled-hex stepping sequence once - // the generator is shared with CPU callers; the benchmark shows the unrolled - // variant winning decisively for divisors >= 131,071. + // Otherwise, find order of 2 mod divisor using the shared unrolled-hex stepping sequence. // TODO: Expose a GPU-first branch here so high divisors leverage the ProcessEightBitWindows // kernel measured fastest in CycleLengthGpuVsCpuBenchmarks, returning the result without // storing it in the shared cache. - ulong order = 1UL; - ulong pow = 2UL; - while (pow != 1UL) - { - pow <<= 1; - if (pow >= divisor) - pow -= divisor; - - order++; - } - - return order; + return divisor.CalculateMersenneDivisorCycleLengthUnrolledHexCpu(); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/ULongExtensions.Cpu.cs b/PerfectNumbers.Core/ULongExtensions.Cpu.cs new file mode 100644 index 0000000..13da576 --- /dev/null +++ b/PerfectNumbers.Core/ULongExtensions.Cpu.cs @@ -0,0 +1,222 @@ +using System; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +public static partial class ULongExtensions +{ + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + ulong modulus = divisor.Modulus; + // This should never happen in production code. If it does, we need to fix it. + // if (exponent == 0UL) + // { + // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + // } + + // There is no point in having extra branch, because in real use-cases it'll barely ever hit. + // Out of all 3_358_944_545 Pow2MontgomeryModWindowedCpu calls, we hit Pow2MontgomeryModSingleBit only 772 times. + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + // } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + + ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); + InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers); + + int index = bitLength - 1; + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MontgomeryMultiply(multiplier, modulus, nPrime); + + index = windowStart - 1; + } + + ThreadStaticPools.UlongPool.Return(oddPowersArray); + + if (!keepMontgomery) + { + result = result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span oddPowers) + { + oddPowers[0] = divisor.MontgomeryTwo; + if (oddPowers.Length == 1) + { + return; + } + + ulong previous; + ulong square = divisor.MontgomeryTwoSquared; + + for (int i = 1; i < oddPowers.Length; i++) + { + previous = oddPowers[i - 1]; + oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) + { + ulong rotationCount = exponent % cycleLength; + return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong CalculateMersenneDivisorCycleLengthUnrolledHexCpu(this ulong divisor) + { + // EvenPerfectBitScanner only routes odd divisors here; keep the guard commented out for benchmarks. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // return 1UL; + // } + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + + if (AdvanceMersenneDivisorCycleStepCpu(ref pow, divisor, ref order)) + { + return order; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool AdvanceMersenneDivisorCycleStepCpu(ref ulong pow, ulong divisor, ref ulong order) + { + pow += pow; + if (pow >= divisor) + { + pow -= divisor; + } + + order++; + return pow == 1UL; + } +} diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs new file mode 100644 index 0000000..895b919 --- /dev/null +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -0,0 +1,131 @@ +using System.Collections.Concurrent; +using System.Runtime.CompilerServices; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core; + +public static partial class ULongExtensions +{ + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulHighGpuCompatible(this ulong x, ulong y) + { + GpuUInt128 product = new(x); + product.Mul64(new GpuUInt128(y)); + return product.High; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) + { + // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, + // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs. + GpuUInt128 state = new(a % modulus); + return state.MulMod(b, modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) + { + // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the + // 5-40× slowdown seen across real-world operand distributions. + GpuUInt128 state = new(a); + return state.MulModWithNativeModulo(b, modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + // We shouldn't hit it in production code. If we do, we're doing something wrong. + // ulong modulus = divisor.Modulus; + // if (exponent == 0UL) + // { + // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + // } + + // We barely ever hit it in production code. It's e.g. 772 calls out of billions + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + // } + + if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult)) + { + return gpuResult; + } + + return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) + { + ulong rotationCount = exponent % cycleLength; + return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false); + } + + private static class Pow2MontgomeryGpuExecutor + { + private static readonly ConcurrentDictionary, MontgomeryDivisorData, byte, ArrayView1D>> KernelCache = new(); + + public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result) + { + result = 0UL; + + GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + + Accelerator accelerator = lease.Accelerator; + // Keep this commented out. It should never happen in production code. + // if (accelerator.AcceleratorType == AcceleratorType.CPU) + // { + // return false; + // } + + var kernel = KernelCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernel); + var launcher = KernelUtil.GetKernel(loaded); + return launcher.CreateLauncherDelegate, MontgomeryDivisorData, byte, ArrayView1D>>(); + }); + + var exponentBuffer = accelerator.Allocate1D(1); + var resultBuffer = accelerator.Allocate1D(1); + + exponentBuffer.View.CopyFromCPU(ref exponent, 1); + // We don't need to worry about any left-overs here. + // resultBuffer.MemSetToZero(); + + // TODO: Have 2 kernels for both cases to have branch-less solution. + byte keepFlag = keepMontgomery ? (byte)1 : (byte)0; + AcceleratorStream stream = lease.Stream; + kernel(stream, 1, exponentBuffer.View, divisor, keepFlag, resultBuffer.View); + stream.Synchronize(); + + resultBuffer.View.CopyToCPU(ref result, 1); + + resultBuffer.Dispose(); + exponentBuffer.Dispose(); + execution.Dispose(); + lease.Dispose(); + // Keep this commented. We don't want to catch any exceptions. All should crash the scanner. + // catch (AcceleratorException) + // { + // return false; + // } + // catch (NotSupportedException) + // { + // return false; + // } + // Intentionally avoid exception handling here; any accelerator failure should crash the scanner. + return true; + } + } +} diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index da55b2d..17b7bea 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -1,14 +1,10 @@ -using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; -using ILGPU; using ILGPU.Algorithms; -using ILGPU.Runtime; -using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; -public static class ULongExtensions +public static partial class ULongExtensions { private const int Pow2WindowSize = 8; private const ulong Pow2WindowFallbackThreshold = 32UL; @@ -459,14 +455,6 @@ public static class ULongExtensions return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulHighGpuCompatible(this ulong x, ulong y) - { - GpuUInt128 product = new(x); - product.Mul64(new GpuUInt128(y)); - return product.High; - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) { @@ -494,24 +482,6 @@ public static class ULongExtensions // MulMod64Benchmarks so CPU callers stop paying for triple modulo operations. public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) - { - // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, - // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs. - GpuUInt128 state = new(a % modulus); - return state.MulMod(b, modulus); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) - { - // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the - // 5-40× slowdown seen across real-world operand distributions. - GpuUInt128 state = new(a); - return state.MulModWithNativeModulo(b, modulus); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) { @@ -528,104 +498,6 @@ public static class ULongExtensions return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) - { - ulong modulus = divisor.Modulus; - // This should never happen in production code. If it does, we need to fix it. - // if (exponent == 0UL) - // { - // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; - // } - - // There is no point in having extra branch, because in real use-cases it'll barely ever hit. - // Out of all 3_358_944_545 Pow2MontgomeryModWindowedCpu calls, we hit Pow2MontgomeryModSingleBit only 772 times. - // if (exponent <= Pow2WindowFallbackThreshold) - // { - // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); - // } - - int bitLength = GetPortableBitLength(exponent); - int windowSize = GetWindowSize(bitLength); - int oddPowerCount = 1 << (windowSize - 1); - - ulong result = divisor.MontgomeryOne; - ulong nPrime = divisor.NPrime; - - ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); - Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); - InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers); - - int index = bitLength - 1; - while (index >= 0) - { - if (((exponent >> index) & 1UL) == 0UL) - { - result = result.MontgomeryMultiply(result, modulus, nPrime); - index--; - continue; - } - - int windowStart = index - windowSize + 1; - if (windowStart < 0) - { - windowStart = 0; - } - - while (((exponent >> windowStart) & 1UL) == 0UL) - { - windowStart++; - } - - int windowLength = index - windowStart + 1; - for (int square = 0; square < windowLength; square++) - { - result = result.MontgomeryMultiply(result, modulus, nPrime); - } - - ulong mask = (1UL << windowLength) - 1UL; - ulong windowValue = (exponent >> windowStart) & mask; - int tableIndex = (int)((windowValue - 1UL) >> 1); - ulong multiplier = oddPowers[tableIndex]; - result = result.MontgomeryMultiply(multiplier, modulus, nPrime); - - index = windowStart - 1; - } - - ThreadStaticPools.UlongPool.Return(oddPowersArray); - - if (!keepMontgomery) - { - result = result.MontgomeryMultiply(1UL, modulus, nPrime); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) - { - // We shouldn't hit it in production code. If we do, we're doing something wrong. - // ulong modulus = divisor.Modulus; - // if (exponent == 0UL) - // { - // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; - // } - - // We barely ever hit it in production code. It's e.g. 772 calls out of billions - // if (exponent <= Pow2WindowFallbackThreshold) - // { - // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); - // } - - if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult)) - { - return gpuResult; - } - - return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static ulong Pow2MontgomeryModWindowedKernel(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { @@ -716,66 +588,6 @@ public static class ULongExtensions return power; } - private static class Pow2MontgomeryGpuExecutor - { - private static readonly ConcurrentDictionary, MontgomeryDivisorData, byte, ArrayView1D>> KernelCache = new(); - - public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result) - { - result = 0UL; - - GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true); - var execution = lease.EnterExecutionScope(); - try - { - Accelerator accelerator = lease.Accelerator; - // Keep this commented out. It should never happen in production code. - // if (accelerator.AcceleratorType == AcceleratorType.CPU) - // { - // return false; - // } - - var kernel = KernelCache.GetOrAdd(accelerator, static accel => - { - var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernel); - var launcher = KernelUtil.GetKernel(loaded); - return launcher.CreateLauncherDelegate, MontgomeryDivisorData, byte, ArrayView1D>>(); - }); - - using var exponentBuffer = accelerator.Allocate1D(1); - using var resultBuffer = accelerator.Allocate1D(1); - - exponentBuffer.View.CopyFromCPU(ref exponent, 1); - // We don't need to worry about any left-overs here. - // resultBuffer.MemSetToZero(); - - // TODO: Have 2 kernels for both cases to have branch-less solution. - byte keepFlag = keepMontgomery ? (byte)1 : (byte)0; - AcceleratorStream stream = lease.Stream; - kernel(stream, 1, exponentBuffer.View, divisor, keepFlag, resultBuffer.View); - stream.Synchronize(); - - resultBuffer.View.CopyToCPU(ref result, 1); - return true; - } - // Keep this commented. We don't want to catch any exceptions. All should crash the scanner. - // catch (AcceleratorException) - // { - // return false; - // } - // catch (NotSupportedException) - // { - // return false; - // } - finally - { - execution.Dispose(); - lease.Dispose(); - } - } - - - } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { @@ -846,51 +658,6 @@ public static class ULongExtensions return Pow2WindowSize; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span oddPowers) - { - oddPowers[0] = divisor.MontgomeryTwo; - if (oddPowers.Length == 1) - { - return; - } - - ulong previous; - ulong square = divisor.MontgomeryTwoSquared; - - for (int i = 1; i < oddPowers.Length; i++) - { - previous = oddPowers[i - 1]; - oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) - { - ulong rotationCount = exponent % cycleLength; - return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) - { - return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) - { - ulong rotationCount = exponent % cycleLength; - return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) - { - return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus) { @@ -1225,4 +992,5 @@ public static class ULongExtensions return false; } + } diff --git a/todo_list.txt b/todo_list.txt index 9f83903..a27c972 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -6,7 +6,6 @@ Divisor Cycle Cache & Generation - EvenPerfectBitScanner/Program.cs:146: Wire GenerateGpu to the unrolled-hex kernel that led the MersenneDivisorCycleLengthGpuBenchmarks once it lands. - EvenPerfectBitScanner/Program.cs:167: Keep a single cached block loaded from disk and honor the configured device when computing ad-hoc cycles for divisors that fall outside that snapshot instead of queuing generation of additional blocks. - EvenPerfectBitScanner/Program.cs:171: Stop reloading the full snapshot once the ad-hoc path streams results straight from the configured device without persisting them, so memory stays bounded while preserving the single cached block strategy. -- PerfectNumbers.Core/MersenneDivisorCycles.cs:209: Replace this naive doubling fallback with the unrolled-hex generator from MersenneDivisorCycleLengthGpuBenchmarks so large divisors avoid the millions of iterations measured in the scalar loop. - PerfectNumbers.Core/MersenneDivisorCycles.cs:212: Route this miss to an ephemeral single-cycle computation on the device selected by the current configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot so we respect the large-p memory limits and avoid extra cache locks. - PerfectNumbers.Core/MersenneDivisorCycles.cs:245: Port this UInt128 path to the unrolled-hex cycle calculator so wide divisors stop relying on the slow shift-and-subtract loop measured in the GPU benchmarks. - PerfectNumbers.Core/MersenneDivisorCycles.cs:464: Reuse the remainder-tracking update from the cycle heuristic so repeated divisions disappear from this loop. @@ -16,7 +15,6 @@ Divisor Cycle Cache & Generation - PerfectNumbers.Core/MersenneDivisorCycles.cs:794: Migrate this generation path to the shared unrolled-hex calculator so the CPU generator matches the GPU benchmark leader for large divisors. - PerfectNumbers.Core/MersenneDivisorCycles.cs:869: Swap these modulo checks for the cached Mod3/Mod5/Mod7/Mod11 helpers so range initialization avoids repeated `%` calls before the unrolled-hex pipeline takes over. - PerfectNumbers.Core/MersenneDivisorCycles.cs:919: Replace this modulo sieve with the cached Mod helpers once the divisor-cycle cache is mandatory so cycle enumeration can drop the slower `%` operations in this hot loop. -- PerfectNumbers.Core/MersenneDivisorCycles.cs:991: Switch this scalar fallback to the unrolled-hex stepping sequence once the generator is shared with CPU callers; the benchmark shows the unrolled variant winning decisively for divisors >= 131,071. - PerfectNumbers.Core/MersenneDivisorCycles.cs:994: Expose a GPU-first branch here so high divisors leverage the ProcessEightBitWindows kernel measured fastest in CycleLengthGpuVsCpuBenchmarks, returning the result without storing it in the shared cache. - PerfectNumbers.Core/UInt128Extensions.cs:372: If the cache lacks this cycle, immediately schedule the configured device (GPU by default) to compute it on the fly and skip inserting it into the cache so wide-order factoring can still leverage cycle stepping without breaching the memory cap or introducing extra synchronization.