From 79695b0c56c1780642d8b9032a0c89e054b9c529 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 14:21:35 +0200 Subject: [PATCH] Switch GPU pow2mod helpers to plain windowed arithmetic --- ...MersenneNumberDivisorByDivisorGpuTester.cs | 38 +-- ...ersenneNumberDivisorResidueGpuEvaluator.cs | 14 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 57 ++-- PerfectNumbers.Core/MersenneDivisorCycles.cs | 9 + PerfectNumbers.Core/ULongExtensions.cs | 261 ++++++++++++++++-- 5 files changed, 305 insertions(+), 74 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index bd96720..ce2c4ec 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -545,10 +545,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { - if (prime <= 1UL) - { - return 0UL; - } + // The by-divisor search emits candidate exponents starting at the first odd prime, so this guard + // would never trigger on the scanner hot path. + // if (prime <= 1UL) + // { + // return 0UL; + // } if (prime - 1UL >= 64UL) { @@ -562,14 +564,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { MontgomeryDivisorData divisor = divisors[index]; ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } + // The divisor staging path filters out non-prime and even moduli before they reach the device, so this guard + // would never fire for EvenPerfectBitScanner workloads. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // hits[index] = 0; + // return; + // } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0; + hits[index] = exponent.Pow2ModWindowedGpu(modulus) == 1UL ? (byte)1 : (byte)0; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -780,15 +784,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } + // Kernel inputs are constrained to odd prime moduli, so this defensive branch would never execute on the configured + // scanning path. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // results[index] = 0UL; + // return; + // } ulong exponent = exponents[index]; - results[index] = exponent.Pow2ModBinaryGpu(modulus); + results[index] = exponent.Pow2ModWindowedGpu(modulus); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index 8e8879a..f14409a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -110,12 +110,14 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable ulong exponent = exponents[idx]; ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[idx] = 0UL; - return; - } + // Residue batches only include valid odd moduli, so the defensive check stays commented to avoid re-validating + // the invariant on every iteration. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // results[idx] = 0UL; + // return; + // } - results[idx] = exponent.Pow2ModBinaryGpu(modulus); + results[idx] = exponent.Pow2ModWindowedGpu(modulus); } } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 2de1d6e..8a572d5 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -805,13 +805,15 @@ internal static partial class PrimeOrderGpuHeuristics for (int i = 0; i < factorCount; i++) { ulong factor = factors[i]; - if (factor <= 1UL) - { - continue; - } + // The factorization routine emits strictly greater-than-one factors for phi on this path, so the + // defensive continue would never execute. + // if (factor <= 1UL) + // { + // continue; + // } ulong reduced = phi / factor; - if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL) + if (reduced.Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { return false; } @@ -827,7 +829,7 @@ internal static partial class PrimeOrderGpuHeuristics if (residue == 1UL || residue == 7UL) { ulong half = phi >> 1; - if (Pow2EqualsOneKernel(half, divisor)) + if ((half).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { order = half; } @@ -861,7 +863,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = order / primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) + if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { order = reduced; continue; @@ -892,7 +894,7 @@ internal static partial class PrimeOrderGpuHeuristics return false; } - if (!Pow2EqualsOneKernel(order, divisor)) + if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL) { return false; } @@ -936,7 +938,7 @@ internal static partial class PrimeOrderGpuHeuristics } reduced /= primeFactor; - if (Pow2EqualsOneKernel(reduced, divisor)) + if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { return false; } @@ -1027,7 +1029,7 @@ internal static partial class PrimeOrderGpuHeuristics } powUsed++; - if (!Pow2EqualsOneKernel(candidate, divisor)) + if ((candidate).Pow2ModWindowedGpu(divisor.Modulus) != 1UL) { continue; } @@ -1327,7 +1329,7 @@ internal static partial class PrimeOrderGpuHeuristics } powUsed++; - if (Pow2EqualsOneKernel(reduced, divisor)) + if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL) { return false; } @@ -1593,11 +1595,6 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) - { - return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL; - } - private static ulong CalculateByDoublingKernel(ulong prime) { ulong order = 1UL; @@ -1698,9 +1695,11 @@ internal static partial class PrimeOrderGpuHeuristics private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) { int length = exponents.Length; + ulong modulus = divisorData.Modulus; + _ = prime; // modulus and divisorData both describe the same q; retain the parameter to preserve the public signature. for (int i = 0; i < length; i++) { - results[i] = Pow2ModCpu(exponents[i], prime, divisorData); + results[i] = exponents[i].Pow2ModWindowedCpu(modulus); } } @@ -1788,16 +1787,6 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorData divisorData) - { - if (modulus <= 1UL) - { - return 0UL; - } - - return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false); - } - private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) { return Pow2ModKernelCache.GetOrAdd(accelerator, static accel => @@ -1811,14 +1800,16 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - remainders[index] = 0UL; - return; - } + // The by-divisor heuristics call this kernel only with odd prime moduli, so the defensive branch stays commented to + // document the invariant without spending a branch in the kernel. + // if (modulus <= 1UL || (modulus & 1UL) == 0UL) + // { + // remainders[index] = 0UL; + // return; + // } ulong exponent = exponents[index]; - remainders[index] = exponent.Pow2ModBinaryGpu(modulus); + remainders[index] = exponent.Pow2ModWindowedGpu(modulus); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 8266171..c918348 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -81,6 +81,8 @@ public class MersenneDivisorCycles if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) { + // Even though production scans never revisit divisors below this bound, the shared cache still serves + // precomputation steps and tests, so keep the fast lookup enabled. var small = Shared._smallCycles; if (small is not null) { @@ -98,6 +100,7 @@ public class MersenneDivisorCycles } Dictionary? cache = factorCache; + // Each by-divisor scan session starts with a null cache and fills it lazily, so keep the allocation guard in place. if (cache is null) { cache = new Dictionary(8); @@ -121,6 +124,8 @@ public class MersenneDivisorCycles if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) { + // The by-divisor path only handles divisors far above this cutoff, but Pollard Rho factoring and unit tests + // still reuse the cached small-cycle table, so keep the fast path enabled. MersenneDivisorCycles shared = Shared; ulong[]? small = shared._smallCycles; if (small is not null) @@ -138,12 +143,16 @@ public class MersenneDivisorCycles private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) { + // The divisor enumeration pipeline has already filtered invalid inputs, but keep the guard so external callers and + // tests continue to receive a defensive answer. if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) { return false; } UInt128 step = (UInt128)exponent << 1; + // Exponents drawn from prime candidates are far below the overflow boundary, yet the broader factoring utilities rely + // on this bail-out when validating arbitrary inputs. if (step == UInt128.Zero) { return false; diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index d8f074c..3cce27d 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -12,6 +12,8 @@ public static class ULongExtensions private const int Pow2WindowSize = 8; private const ulong Pow2WindowFallbackThreshold = 32UL; private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); + private const int Pow2WindowSizeGpu = 5; + private const int Pow2WindowOddPowerCountGpu = 1 << (Pow2WindowSizeGpu - 1); public static ulong CalculateOrder(this ulong q) { @@ -310,25 +312,26 @@ public static class ULongExtensions return state.MulModWithNativeModulo(b, modulus); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2ModBinaryCpu(this ulong exponent, ulong modulus) - { - return Pow2ModWindowedCpu(exponent, modulus); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus) { - if (modulus <= 1UL) - { - return 0UL; - } + // The by-divisor pipeline only schedules odd prime moduli, so the defensive check below + // would never trigger on that path. Leaving it commented documents the invariant without + // forcing a branch in the hot loop. + // if (modulus <= 1UL) + // { + // return 0UL; + // } - if (exponent == 0UL) - { - return 1UL % modulus; - } + // Exponents flow from the strictly positive candidate list on the by-divisor scanner, so + // this guard would be redundant there. + // if (exponent == 0UL) + // { + // return 1UL % modulus; + // } + // Even though the by-divisor pipeline only feeds very large exponents, the legacy callers and unit tests + // still hit this small-exponent fast path, so keep the fallback in place. if (exponent <= Pow2WindowFallbackThreshold) { return Pow2ModBinaryFallback(exponent, modulus); @@ -341,11 +344,12 @@ public static class ULongExtensions Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : new ulong[oddPowerCount]; - // Reuse oddPowerStorage to expose only the populated slice and avoid another local span. oddPowerStorage = oddPowerStorage[..oddPowerCount]; InitializePlainOddPowers(modulus, oddPowerStorage); - ulong result = 1UL % modulus; + // With modulus guaranteed to exceed one on the scanning path, the modulo operation would + // always yield one, so keep the literal to highlight the intended base value. + ulong result = 1UL; int index = bitLength - 1; while (index >= 0) @@ -387,6 +391,8 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong Pow2ModBinaryFallback(ulong exponent, ulong modulus) { + // The scanning path keeps modulus strictly greater than one, but other utility callers still + // rely on this fallback for tiny exponents, so leave the modulo in place for safety. ulong remainingExponent = exponent; ulong result = 1UL % modulus; ulong baseValue = 2UL % modulus; @@ -411,13 +417,77 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus) + public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) { - if (modulus <= 1UL) + // The GPU pipeline feeds only odd prime moduli, so the removed guard would never fire for + // EvenPerfectBitScanner workloads. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + + // Exponents fed into the GPU kernels are positive on the scanning path, yet the zero case remains for + // completeness when unit tests probe boundary behavior. + if (exponent == 0UL) { - return 0UL; + return 1UL; } + // The GPU kernels mostly process large exponents as well, but the binary fallback is still useful when + // unit tests probe edge cases around the threshold. + if (exponent <= Pow2WindowFallbackThreshold) + { + return Pow2ModBinaryGpuFallback(exponent, modulus); + } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetGpuWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + PlainOddPowerTableGpu oddPowers = new PlainOddPowerTableGpu(modulus, oddPowerCount); + + ulong result = 1UL; + int index = bitLength - 1; + + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MulMod64GpuCompatible(result, modulus); + 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.MulMod64GpuCompatible(result, modulus); + } + + ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL); + int tableIndex = (int)((windowValue - 1UL) >> 1); + result = result.MulMod64GpuCompatible(oddPowers.Get(tableIndex), modulus); + + index = windowStart - 1; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2ModBinaryGpuFallback(ulong exponent, ulong modulus) + { ulong remainingExponent = exponent; ulong result = 1UL % modulus; ulong baseValue = 2UL % modulus; @@ -678,10 +748,19 @@ public static class ULongExtensions return Pow2WindowSize; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetGpuWindowSize(int bitLength) + { + int window = GetWindowSize(bitLength); + return window > Pow2WindowSizeGpu ? Pow2WindowSizeGpu : window; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void InitializePlainOddPowers(ulong modulus, Span destination) { int oddPowerCount = destination.Length; + // The by-divisor kernels always request at least one odd power, but other helpers reuse this initializer + // for smaller workloads, so the early exit stays in place. if (oddPowerCount == 0) { return; @@ -703,6 +782,150 @@ public static class ULongExtensions } } + private struct PlainOddPowerTableGpu + { + private readonly int _count; + private ulong _value0; + private ulong _value1; + private ulong _value2; + private ulong _value3; + private ulong _value4; + private ulong _value5; + private ulong _value6; + private ulong _value7; + private ulong _value8; + private ulong _value9; + private ulong _value10; + private ulong _value11; + private ulong _value12; + private ulong _value13; + private ulong _value14; + private ulong _value15; + + public PlainOddPowerTableGpu(ulong modulus, int requestedCount) + { + _count = requestedCount > Pow2WindowOddPowerCountGpu ? Pow2WindowOddPowerCountGpu : requestedCount; + _value0 = 0UL; + _value1 = 0UL; + _value2 = 0UL; + _value3 = 0UL; + _value4 = 0UL; + _value5 = 0UL; + _value6 = 0UL; + _value7 = 0UL; + _value8 = 0UL; + _value9 = 0UL; + _value10 = 0UL; + _value11 = 0UL; + _value12 = 0UL; + _value13 = 0UL; + _value14 = 0UL; + _value15 = 0UL; + + if (_count == 0) + { + return; + } + + ulong baseValue = 2UL % modulus; + _value0 = baseValue; + if (_count == 1) + { + return; + } + + ulong square = baseValue.MulMod64GpuCompatible(baseValue, modulus); + ulong current = baseValue; + for (int i = 1; i < _count; i++) + { + current = current.MulMod64GpuCompatible(square, modulus); + SetValue(i, current); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong Get(int index) + { + if ((uint)index >= (uint)_count) + { + return 0UL; + } + + return index switch + { + 0 => _value0, + 1 => _value1, + 2 => _value2, + 3 => _value3, + 4 => _value4, + 5 => _value5, + 6 => _value6, + 7 => _value7, + 8 => _value8, + 9 => _value9, + 10 => _value10, + 11 => _value11, + 12 => _value12, + 13 => _value13, + 14 => _value14, + _ => _value15, + }; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void SetValue(int index, ulong value) + { + switch (index) + { + case 1: + _value1 = value; + break; + case 2: + _value2 = value; + break; + case 3: + _value3 = value; + break; + case 4: + _value4 = value; + break; + case 5: + _value5 = value; + break; + case 6: + _value6 = value; + break; + case 7: + _value7 = value; + break; + case 8: + _value8 = value; + break; + case 9: + _value9 = value; + break; + case 10: + _value10 = value; + break; + case 11: + _value11 = value; + break; + case 12: + _value12 = value; + break; + case 13: + _value13 = value; + break; + case 14: + _value14 = value; + break; + case 15: + _value15 = value; + break; + } + } + } + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span destination) { int oddPowerCount = destination.Length;