From 79e8706644b83c64044a5ecee108f883dc866572 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 19 Oct 2025 20:21:09 +0200 Subject: [PATCH] Restore original guard comments in pow2 helpers --- .../MulMod64Benchmarks.cs | 8 +- .../Gpu/ULongExtensionsGpuTests.cs | 64 +++++ .../Gpu/Kernels/DivisorByDivisorKernels.cs | 5 +- .../Gpu/Kernels/Pow2MontgomeryKernels.cs | 3 +- .../Gpu/Kernels/PrimeOrderKernels.cs | 7 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 3 +- PerfectNumbers.Core/ULongExtensions.Cpu.cs | 51 +++- PerfectNumbers.Core/ULongExtensions.Gpu.cs | 249 +++++++++++++++++- PerfectNumbers.Core/ULongExtensions.cs | 135 +--------- 9 files changed, 367 insertions(+), 158 deletions(-) create mode 100644 PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index 388cdcd..b995d0e 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -157,7 +157,7 @@ public class MulMod64Benchmarks } /// - /// Uses the GPU-friendly helper; dense inputs now + /// Uses the GPU-friendly helper; dense inputs now /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at /// 3.59 ns. /// @@ -170,11 +170,11 @@ public class MulMod64Benchmarks { // TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity // is not required; that baseline stayed 6–17× faster for dense 64-bit operands. - return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus); + return Input.Left.MulMod64Gpu(Input.Right, Input.Modulus); } /// - /// Uses the deferred native-modulo helper ; + /// Uses the deferred native-modulo helper ; /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns. /// /// @@ -186,7 +186,7 @@ public class MulMod64Benchmarks { // TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks; // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes. - return Input.Left.MulMod64GpuCompatibleDeferred(Input.Right, Input.Modulus); + return Input.Left.MulMod64GpuDeferred(Input.Right, Input.Modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs new file mode 100644 index 0000000..0f9f674 --- /dev/null +++ b/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs @@ -0,0 +1,64 @@ +using FluentAssertions; +using PerfectNumbers.Core; +using Xunit; + +namespace PerfectNumbers.Core.Tests.Gpu; + +public class ULongExtensionsGpuTests +{ + [Theory] + [Trait("Category", "Fast")] + [InlineData(1UL, 3UL)] + [InlineData(2UL, 5UL)] + [InlineData(3UL, 7UL)] + [InlineData(63UL, 97UL)] + [InlineData(512UL, 12289UL)] + [InlineData(1024UL, 196613UL)] + [InlineData(4096UL, 999983UL)] + [InlineData(65535UL, 2147483647UL)] + public void Pow2ModWindowedGpu_matches_montgomery_cpu_for_known_values(ulong exponent, ulong modulus) + { + MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + ulong expected = exponent.Pow2MontgomeryModWindowedCpu(divisor, keepMontgomery: false); + + exponent.Pow2ModWindowedGpu(modulus).Should().Be(expected); + } + + [Fact] + [Trait("Category", "Fast")] + public void Pow2ModWindowedGpu_matches_cpu_for_sampled_inputs() + { + ulong[] moduli = + { + 3UL, + 5UL, + 97UL, + 12289UL, + 196613UL, + 4294967291UL, + 9223372036854775783UL, + }; + + ulong[] exponents = + { + 1UL, + 2UL, + 7UL, + 32UL, + 63UL, + 2048UL, + 65535UL, + }; + + foreach (ulong modulus in moduli) + { + MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + foreach (ulong exponent in exponents) + { + ulong expected = exponent.Pow2MontgomeryModWindowedCpu(divisor, keepMontgomery: false); + exponent.Pow2ModWindowedGpu(modulus).Should().Be(expected); + } + } + } + +} diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index f3c0a17..eb3d977 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -1,5 +1,6 @@ using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -16,7 +17,7 @@ internal static class DivisorByDivisorKernels } ulong exponent = exponents[index]; - hits[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; + hits[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponent, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; } public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) @@ -148,6 +149,6 @@ internal static class DivisorByDivisorKernels ulong exponent = exponents[index]; - results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true); + results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponent, keepMontgomery: true); } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs index 806e828..64120f7 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs @@ -1,5 +1,6 @@ using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -12,6 +13,6 @@ internal static class Pow2MontgomeryKernels byte keepMontgomery, ArrayView1D results) { - results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0); + results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponents[index], keepMontgomery != 0); } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs index c7076f5..988e603 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs @@ -4,6 +4,7 @@ using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Algorithms.Random; using ILGPU.Runtime; +using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; @@ -346,7 +347,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = phi / factor; - if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL) + if (ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, reduced, keepMontgomery: false) == 1UL) { return false; } @@ -1130,7 +1131,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) { - return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL; + return ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponent, keepMontgomery: false) == 1UL; } private static ulong CalculateByDoublingKernel(ulong prime) @@ -1155,7 +1156,7 @@ internal static partial class PrimeOrderGpuHeuristics internal static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { ulong exponent = exponents[index]; - remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false); + remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponent, keepMontgomery: false); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index a855dd7..716221a 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -3,6 +3,7 @@ using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; using System.Numerics; +using PerfectNumbers.Core; using System.Runtime.InteropServices; using System.Runtime.CompilerServices; using ILGPU.Runtime.OpenCL; @@ -651,7 +652,7 @@ internal static partial class PrimeOrderGpuHeuristics return 0UL; } - return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisorData, keepMontgomery: false); + return ULongExtensions.Pow2MontgomeryModWindowedGpu(divisorData, exponent, keepMontgomery: false); } private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core/ULongExtensions.Cpu.cs b/PerfectNumbers.Core/ULongExtensions.Cpu.cs index 13da576..3de0b55 100644 --- a/PerfectNumbers.Core/ULongExtensions.Cpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Cpu.cs @@ -1,4 +1,5 @@ using System; +using System.Numerics; using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; @@ -9,7 +10,7 @@ public static partial class ULongExtensions 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. + // This should never happen in production code. // if (exponent == 0UL) // { // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; @@ -22,8 +23,8 @@ public static partial class ULongExtensions // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); // } - int bitLength = GetPortableBitLength(exponent); - int windowSize = GetWindowSize(bitLength); + int bitLength = GetPortableBitLengthCpu(exponent); + int windowSize = GetWindowSizeCpu(bitLength); int oddPowerCount = 1 << (windowSize - 1); ulong result = divisor.MontgomeryOne; @@ -219,4 +220,48 @@ public static partial class ULongExtensions order++; return pow == 1UL; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetPortableBitLengthCpu(ulong value) + { + // Keep this commented out. It will never happen in production code. + // if (value == 0UL) + // { + // return 0; + // } + + return 64 - BitOperations.LeadingZeroCount(value); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetWindowSizeCpu(int bitLength) + { + if (bitLength <= 6) + { + 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 Pow2WindowSize; + } + } diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index 895b919..cd445d1 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -1,6 +1,7 @@ using System.Collections.Concurrent; using System.Runtime.CompilerServices; using ILGPU; +using ILGPU.Algorithms; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; @@ -17,7 +18,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) + public static ulong MulMod64Gpu(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. @@ -26,7 +27,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) + public static ulong MulMod64GpuDeferred(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. @@ -34,6 +35,147 @@ public static partial class ULongExtensions return state.MulModWithNativeModulo(b, modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) + { + // EvenPerfectBitScanner never routes moduli ≤ 1 through this path; keep the guard commented out so the hot loop stays branch-free. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + + // EvenPerfectBitScanner always works with positive exponents here; leave the defensive check commented for future reference. + // if (exponent == 0UL) + // { + // return 1UL % modulus; + // } + + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + // oddPowerCount is always ≥ 1 for the production workloads handled by EvenPerfectBitScanner. + // if (oddPowerCount <= 0) + // { + // return 1UL; + // } + + ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); + + InitializeStandardOddPowers(modulus, oddPowers); + + ulong result = 1UL; + int index = bitLength - 1; + + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + ulong squared = result.MulMod64Gpu(result, modulus); + result = currentBit == 0UL ? squared : result; + index = currentBit == 0UL ? index - 1 : index; + if (currentBit == 0UL) + { + continue; + } + + int windowStart = index - windowSize + 1; + // EvenPerfectBitScanner guarantees windowStart stays non-negative for production exponents, so keep the guard commented out to avoid extra branching. + // if (windowStart < 0) + // { + // windowStart = 0; + // } + + windowStart = GetNextSetBitIndexGpu(exponent, windowStart); + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MulMod64Gpu(result, modulus); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MulMod64Gpu(multiplier, modulus); + + index = windowStart - 1; + } + + ThreadStaticPools.UlongPool.Return(oddPowersArray); + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong Pow2MontgomeryModWindowedGpu(in MontgomeryDivisorData divisor, ulong exponent, bool keepMontgomery) + { + ulong modulus = divisor.Modulus; + // Reuse the optimized MontgomeryMultiply extension highlighted in MontgomeryMultiplyBenchmarks so the GPU path matches the fastest CPU reduction. + // We shouldn't hit it in production code. If we do, we're doing something wrong. + // 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); + // } + + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + + int index = bitLength - 1; + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + ulong squared = MontgomeryMultiply(result, result, modulus, nPrime); + result = currentBit == 0UL ? squared : result; + index -= (int)(currentBit ^ 1UL); + if (currentBit == 0UL) + { + continue; + } + + int windowStart = index - windowSize + 1; + // if (windowStart < 0) + // { + // windowStart = 0; + // } + // Clamp the negative offset without branching so the GPU loop stays divergence-free. + // This still handles windows that would otherwise extend past the most significant bit. + int negativeMask = windowStart >> 31; + windowStart &= ~negativeMask; + + windowStart = GetNextSetBitIndexGpu(exponent, windowStart); + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = MontgomeryMultiply(result, result, modulus, nPrime); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + ulong multiplier = ComputeMontgomeryOddPowerGpu(windowValue, divisor, modulus, nPrime); + result = MontgomeryMultiply(result, multiplier, modulus, nPrime); + + index = windowStart - 1; + } + + if (!keepMontgomery) + { + result = MontgomeryMultiply(result, 1UL, modulus, nPrime); + } + + return result; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { @@ -50,12 +192,97 @@ public static partial class ULongExtensions // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); // } - if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult)) + return Pow2MontgomeryGpuExecutor.Execute(exponent, divisor, keepMontgomery); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializeStandardOddPowers(ulong modulus, Span oddPowers) + { + // oddPowers spans at least one entry for production workloads; keep the guard commented out to avoid redundant checks in the hot loop. + // if (oddPowers.IsEmpty) + // { + // return; + // } + + // EvenPerfectBitScanner feeds odd prime moduli (≥ 3) here, so the base value stays within range without a modulo reduction. + ulong baseValue = 2UL; + oddPowers[0] = baseValue; + if (oddPowers.Length == 1) { - return gpuResult; + return; } - return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery); + ulong square = baseValue.MulMod64Gpu(baseValue, modulus); + for (int i = 1; i < oddPowers.Length; i++) + { + ulong previous = oddPowers[i - 1]; + oddPowers[i] = previous.MulMod64Gpu(square, modulus); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeMontgomeryOddPowerGpu(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime) + { + ulong baseValue = divisor.MontgomeryTwo; + ulong power = divisor.MontgomeryOne; + ulong remaining = exponent; + + while (remaining != 0UL) + { + ulong bit = remaining & 1UL; + ulong multiplied = MontgomeryMultiply(power, baseValue, modulus, nPrime); + ulong mask = (ulong)-(long)bit; + power = (power & ~mask) | (multiplied & mask); + + remaining >>= 1; + if (remaining == 0UL) + { + break; + } + + baseValue = MontgomeryMultiply(baseValue, baseValue, modulus, nPrime); + } + + return power; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetPortableBitLengthGpu(ulong value) + { + // Keep this commented out. It will never happen in production code. + // if (value == 0UL) + // { + // return 0; + // } + + return 64 - XMath.LeadingZeroCount(value); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetWindowSizeGpu(int n) + { + int w = Pow2WindowSize; + + w = (n <= 671) ? 7 : w; + w = (n <= 239) ? 6 : w; + w = (n <= 79) ? 5 : w; + w = (n <= 23) ? 4 : w; + + int m = (n >= 1) ? n : 1; + w = (n <= 6) ? m : w; + + return w; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetNextSetBitIndexGpu(ulong exponent, int startIndex) + { + ulong guard = (ulong)(((long)startIndex - 64) >> 63); + int shift = startIndex & 63; + ulong mask = (~0UL << shift) & guard; + ulong masked = exponent & mask; + + return XMath.TrailingZeroCount(masked); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -75,9 +302,9 @@ public static partial class ULongExtensions { private static readonly ConcurrentDictionary, MontgomeryDivisorData, byte, ArrayView1D>> KernelCache = new(); - public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result) + public static ulong Execute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { - result = 0UL; + ulong result = 0UL; GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true); var execution = lease.EnterExecutionScope(); @@ -86,7 +313,7 @@ public static partial class ULongExtensions // Keep this commented out. It should never happen in production code. // if (accelerator.AcceleratorType == AcceleratorType.CPU) // { - // return false; + // return 0UL; // } var kernel = KernelCache.GetOrAdd(accelerator, static accel => @@ -118,14 +345,14 @@ public static partial class ULongExtensions // Keep this commented. We don't want to catch any exceptions. All should crash the scanner. // catch (AcceleratorException) // { - // return false; + // return 0UL; // } // catch (NotSupportedException) // { - // return false; + // return 0UL; // } // Intentionally avoid exception handling here; any accelerator failure should crash the scanner. - return true; + return result; } } } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 17b7bea..33c7585 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -498,98 +498,8 @@ public static partial class ULongExtensions return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2MontgomeryModWindowedKernel(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) - { - ulong modulus = divisor.Modulus; - // This should never happen in production code - // if (exponent == 0UL) - // { - // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; - // } - - // We barely every hit it from production code. It's only 772 calls out of billions. - // if (exponent <= Pow2WindowFallbackThreshold) - // { - // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); - // } - - int bitLength = GetPortableBitLength(exponent); - int windowSize = GetWindowSize(bitLength); - ulong result = divisor.MontgomeryOne; - ulong nPrime = divisor.NPrime; - - 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; - ulong multiplier = ComputeMontgomeryOddPower(windowValue, divisor, modulus, nPrime); - result = result.MontgomeryMultiply(multiplier, modulus, nPrime); - - index = windowStart - 1; - } - - if (!keepMontgomery) - { - result = result.MontgomeryMultiply(1UL, modulus, nPrime); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeMontgomeryOddPower(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime) - { - ulong baseValue = divisor.MontgomeryTwo; - ulong power = divisor.MontgomeryOne; - ulong remaining = exponent; - - while (remaining != 0UL) - { - if ((remaining & 1UL) != 0UL) - { - power = power.MontgomeryMultiply(baseValue, modulus, nPrime); - } - - remaining >>= 1; - if (remaining == 0UL) - { - break; - } - - baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); - } - - return power; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { ulong modulus = divisor.Modulus; ulong nPrime = divisor.NPrime; @@ -616,47 +526,6 @@ public static partial class ULongExtensions return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetPortableBitLength(ulong value) - { - // Keep this commented out. It will never happen in production code. - // if (value == 0UL) - // { - // return 0; - // } - - return 64 - XMath.LeadingZeroCount(value); - } - - private static int GetWindowSize(int bitLength) - { - if (bitLength <= 6) - { - 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 Pow2WindowSize; - } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus)