From 7d112c347f9a20eb9a4d615e9c421be44f9a96e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 22 Sep 2025 23:24:19 +0200 Subject: [PATCH] Accelerate Mersenne modular multiplication paths --- EvenPerfectBitScanner/Program.cs | 4 +- PerfectNumbers.Core/Gpu/GpuUInt128.cs | 119 +++++++++++++++- PerfectNumbers.Core/UInt128Extensions.cs | 12 +- PerfectNumbers.Core/UIntExtensions.cs | 87 +----------- PerfectNumbers.Core/ULongExtensions.cs | 173 +++++++++-------------- 5 files changed, 196 insertions(+), 199 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 67d9a7a..ff2676a 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -520,7 +520,7 @@ internal static class Program RleBlacklist.Load(_rleBlacklistPath!); } - ulong remainder = currentP.Mod6(); + ulong remainder = currentP % 6UL; if (currentP == InitialP && string.IsNullOrEmpty(filterFile)) { // bool passedAllTests = IsEvenPerfectCandidate(InitialP, out bool searchedMersenne, out bool detailedCheck); @@ -1649,7 +1649,7 @@ internal static class Program // Optional: RLE blacklist and binary-threshold filters on p (safe only when configured) if (p <= _rleHardMaxP) { - if (!_rleOnlyLast7 || p.Mod10() == 7UL) + if (!_rleOnlyLast7 || p % 10UL == 7UL) { if (RleBlacklist.IsLoaded() && RleBlacklist.Matches(p)) { diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index 725780f..1b1487c 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -1057,13 +1057,27 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong MulMod64(ulong a, ulong b, ulong modulus) { - // Note: 64-bit % is supported in ILGPU kernels, but the straightforward - // (UInt128)a*b % modulus is not, because 128-bit multiply is not - // supported in device code. We therefore implement a branch-free - // shift-add reduction here based solely on 64-bit ops. + if (modulus <= 1UL) + { + return 0UL; + } + + ulong aReduced = a % modulus; + ulong bReduced = b % modulus; + if (aReduced == 0UL || bReduced == 0UL) + { + return 0UL; + } + + if (TryGetMersenneExponent(modulus, out int exponent)) + { + return MulModMersenne(aReduced, bReduced, modulus, exponent); + } + + // Fallback for general moduli using the branch-free shift-add reducer. ulong result = 0UL; - ulong x = a % modulus; - ulong y = b; + ulong x = aReduced; + ulong y = bReduced; while (y != 0UL) { if ((y & 1UL) != 0UL) @@ -1087,6 +1101,99 @@ public struct GpuUInt128 : IComparable, IEquatable return result; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong MulModMersenne(ulong a, ulong b, ulong modulus, int exponent) + { + Mul64Parts(a, b, out ulong productHigh, out ulong productLow); + ulong mask = exponent == 64 ? ulong.MaxValue : modulus; + + ulong currentHigh = productHigh; + ulong currentLow = productLow; + + do + { + ulong lower = currentLow & mask; + ShiftRight128(currentHigh, currentLow, exponent, out ulong shiftedHigh, out ulong shiftedLow); + currentLow = lower + shiftedLow; + currentHigh = shiftedHigh; + if (currentLow < lower) + { + currentHigh++; + } + } + while (currentHigh != 0UL); + + ulong result = currentLow; + while (result >= modulus) + { + result -= modulus; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool TryGetMersenneExponent(ulong modulus, out int exponent) + { + if (modulus <= 1UL) + { + exponent = 0; + return false; + } + + ulong plusOne = modulus + 1UL; + if (plusOne == 0UL) + { + exponent = 64; + return true; + } + + if ((plusOne & (plusOne - 1UL)) != 0UL) + { + exponent = 0; + return false; + } + + int bits = 0; + while ((plusOne & 1UL) == 0UL) + { + plusOne >>= 1; + bits++; + } + + exponent = bits; + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void ShiftRight128(ulong high, ulong low, int shift, out ulong newHigh, out ulong newLow) + { + if (shift == 0) + { + newHigh = high; + newLow = low; + return; + } + + if (shift < 64) + { + newLow = (low >> shift) | (high << (64 - shift)); + newHigh = high >> shift; + return; + } + + if (shift == 64) + { + newLow = high; + newHigh = 0UL; + return; + } + + int extra = shift - 64; + newLow = high >> extra; + newHigh = 0UL; + } + // Montgomery core for 64-bit operands. Returns a*b*R^{-1} mod modulus. [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong MontMul64(ulong aR, ulong bR, ulong modulus, ulong nPrime) diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 99fda8a..ee50019 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -218,11 +218,11 @@ public static class UInt128Extensions while (high != zero) { // 2^64 ≡ 6 (mod 10) - result = (result + (ulong)high * 6UL).Mod10(); + result = (result + (ulong)high * 6UL) % 10UL; high >>= 64; } - return result.Mod10(); + return result % 10UL; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -231,14 +231,14 @@ public static class UInt128Extensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod3(this UInt128 value) { - ulong remainder = ((ulong)value).Mod3() + ((ulong)(value >> 64)).Mod3(); + ulong remainder = ((ulong)value % 3UL) + ((ulong)(value >> 64) % 3UL); return remainder >= 3UL ? remainder - 3UL : remainder; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod5(this UInt128 value) { - ulong remainder = ((ulong)value).Mod5() + ((ulong)(value >> 64)).Mod5(); + ulong remainder = ((ulong)value % 5UL) + ((ulong)(value >> 64) % 5UL); return remainder >= 5UL ? remainder - 5UL : remainder; } @@ -259,7 +259,7 @@ public static class UInt128Extensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod7(this UInt128 value) { - ulong remainder = ((ulong)value).Mod7() + ((ulong)(value >> 64)).Mod7() * 2UL; + ulong remainder = ((ulong)value % 7UL) + ((ulong)(value >> 64) % 7UL) * 2UL; while (remainder >= 7UL) { remainder -= 7UL; @@ -271,7 +271,7 @@ public static class UInt128Extensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod11(this UInt128 value) { - ulong remainder = ((ulong)value).Mod11() + ((ulong)(value >> 64)).Mod11() * 5UL; + ulong remainder = ((ulong)value % 11UL) + ((ulong)(value >> 64) % 11UL) * 5UL; while (remainder >= 11UL) { remainder -= 11UL; diff --git a/PerfectNumbers.Core/UIntExtensions.cs b/PerfectNumbers.Core/UIntExtensions.cs index 3706e18..ba23e63 100644 --- a/PerfectNumbers.Core/UIntExtensions.cs +++ b/PerfectNumbers.Core/UIntExtensions.cs @@ -4,99 +4,26 @@ namespace PerfectNumbers.Core; public static class UIntExtensions { - private static readonly byte[] Mod6Lookup = { 0, 3, 4, 1, 2, 5 }; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static uint Mod3(this uint value) => value % 3U; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod3(this uint value) - { - ulong wideValue = 0UL; - uint quotient = 0U; - uint remainder = 0U; - - wideValue = (ulong)value * 0xAAAAAAABUL; - quotient = (uint)(wideValue >> 33); - remainder = value - quotient * 3U; - - return remainder; - } + public static uint Mod5(this uint value) => value % 5U; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod5(this uint value) - { - ulong wideValue = 0UL; - uint quotient = 0U; - uint remainder = 0U; - - wideValue = (ulong)value * 0xCCCCCCCDUL; - quotient = (uint)(wideValue >> 34); - remainder = value - quotient * 5U; - - return remainder; - } + public static uint Mod6(this uint value) => value % 6U; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod6(this uint value) - { - return Mod6Lookup[(int)(((value.Mod3() << 1) | (value & 1U)))]; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod7(this uint value) - { - uint remainder = 0U; - uint temp = 0U; - uint chunk = 0U; - - temp = value; - while (temp != 0U) - { - chunk = temp & 7U; - remainder += chunk; - temp >>= 3; - if (remainder >= 7U) - { - remainder -= 7U; - } - } - - return remainder >= 7U ? remainder - 7U : remainder; - } + public static uint Mod7(this uint value) => value % 7U; [MethodImpl(MethodImplOptions.AggressiveInlining)] public static uint Mod8(this uint value) => value & 7U; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod10(this uint value) - { - ulong wideValue = 0UL; - uint quotient = 0U; - uint remainder = 0U; - - wideValue = (ulong)value * 0xCCCCCCCDUL; - quotient = (uint)(wideValue >> 35); - remainder = value - quotient * 10U; - - return remainder; - } + public static uint Mod10(this uint value) => value % 10U; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod11(this uint value) - { - uint remainder = 0U; - uint temp = 0U; - uint chunk = 0U; - - temp = value; - while (temp != 0U) - { - chunk = temp & 1023U; - remainder += chunk; - temp >>= 10; - remainder -= 11U * (remainder / 11U); - } - - return remainder; - } + public static uint Mod11(this uint value) => value % 11U; [MethodImpl(MethodImplOptions.AggressiveInlining)] public static uint Mod128(this uint value) => value & 127U; diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 61fa6f4..a311acb 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -57,7 +57,6 @@ public static class ULongExtensions return q; } - private static readonly byte[] Mod6Lookup = { 0, 3, 4, 1, 2, 5 }; private static readonly byte[] Mod7ByteCoefficients = { 1, 4, 2, 1, 4, 2, 1, 4 }; private static readonly byte[] Mod11ByteCoefficients = { 1, 3, 9, 5, 4, 1, 3, 9 }; @@ -187,121 +186,25 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod10(this ulong value) - { - ulong mod5 = 0UL; - ulong result = 0UL; - - mod5 = value.Mod5(); - - if ((value & 1UL) == 0UL) - { - result = mod5 switch - { - 0UL => 0UL, - 1UL => 6UL, - 2UL => 2UL, - 3UL => 8UL, - _ => 4UL, - }; - } - else - { - result = mod5 switch - { - 0UL => 5UL, - 1UL => 1UL, - 2UL => 7UL, - 3UL => 3UL, - _ => 9UL, - }; - } - - return result; - } + public static ulong Mod10(this ulong value) => value % 10UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod8(this ulong value) => value & 7UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod3(this ulong value) - { - uint sum = 0U; - ulong temp = 0UL; - - temp = value & 0xFFFFUL; - sum = (uint)temp; - temp = (value >> 16) & 0xFFFFUL; - sum += (uint)temp; - temp = (value >> 32) & 0xFFFFUL; - sum += (uint)temp; - temp = (value >> 48) & 0xFFFFUL; - sum += (uint)temp; - - return sum.Mod3(); - } + public static ulong Mod3(this ulong value) => value % 3UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod5(this ulong value) - { - uint sum = 0U; - ulong temp = 0UL; - - temp = value & 0xFFFFUL; - sum = (uint)temp; - temp = (value >> 16) & 0xFFFFUL; - sum += (uint)temp; - temp = (value >> 32) & 0xFFFFUL; - sum += (uint)temp; - temp = (value >> 48) & 0xFFFFUL; - sum += (uint)temp; - - return sum.Mod5(); - } + public static ulong Mod5(this ulong value) => value % 5UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod6(this ulong value) - { - return Mod6Lookup[(int)(((value.Mod3() << 1) | (value & 1UL)))]; - } + public static ulong Mod6(this ulong value) => value % 6UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod7(this ulong value) - { - uint low = 0U; - uint high = 0U; - ulong remainder = 0UL; - - low = (uint)value; - high = (uint)(value >> 32); - - remainder = low.Mod7() + (ulong)high.Mod7() * 4UL; - while (remainder >= 7UL) - { - remainder -= 7UL; - } - - return remainder; - } + public static ulong Mod7(this ulong value) => value % 7UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod11(this ulong value) - { - uint low = 0U; - uint high = 0U; - ulong remainder = 0UL; - - low = (uint)value; - high = (uint)(value >> 32); - - remainder = low.Mod11() + (ulong)high.Mod11() * 4UL; - while (remainder >= 11UL) - { - remainder -= 11UL; - } - - return remainder; - } + public static ulong Mod11(this ulong value) => value % 11UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod128(this ulong value) => value & 127UL; @@ -377,8 +280,68 @@ public static class ULongExtensions return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64(this ulong a, ulong b, ulong modulus) + { + if (modulus <= 1UL) + { + return 0UL; + } + + if (TryGetMersenneExponent(modulus, out int exponent)) + { + return MulModMersenne(a, b, modulus, exponent); + } + + return (ulong)(((UInt128)a * b) % modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool TryGetMersenneExponent(ulong modulus, out int exponent) + { + if (modulus == 0UL) + { + exponent = 0; + return false; + } + + ulong plusOne = modulus + 1UL; + if (plusOne == 0UL) + { + exponent = 64; + return true; + } + + if ((plusOne & (plusOne - 1UL)) != 0UL) + { + exponent = 0; + return false; + } + + exponent = BitOperations.TrailingZeroCount(plusOne); + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong MulModMersenne(ulong a, ulong b, ulong modulus, int exponent) + { + UInt128 mask = (UInt128)modulus; + UInt128 product = (UInt128)a * b; + UInt128 folded = (product & mask) + (product >> exponent); + + while (folded > mask) + { + folded = (folded & mask) + (folded >> exponent); + } + + ulong result = (ulong)folded; + if (result >= modulus) + { + result -= modulus; + } + + return result; + } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus)