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)