Restore original guard comments in pow2 helpers
This commit is contained in:
1 parent
da7a6cb23d
commit
79e8706644
9 files changed
+366
-157
No files matched your search
@@ -157,7 +157,7 @@ public class MulMod64Benchmarks
|
|||||||
}
|
}
|
||||||
|
|
||||||
/// <summary>
|
/// <summary>
|
||||||
/// Uses the GPU-friendly <see cref="ULongExtensions.MulMod64GpuCompatible(ulong, ulong, ulong)"/> helper; dense inputs now
|
/// Uses the GPU-friendly <see cref="ULongExtensions.MulMod64Gpu(ulong, ulong, ulong)"/> 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
|
/// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at
|
||||||
/// 3.59 ns.
|
/// 3.59 ns.
|
||||||
/// </summary>
|
/// </summary>
|
||||||
@@ -170,11 +170,11 @@ public class MulMod64Benchmarks
|
|||||||
{
|
{
|
||||||
// TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity
|
// 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.
|
// 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);
|
||||||
}
|
}
|
||||||
|
|
||||||
/// <summary>
|
/// <summary>
|
||||||
/// Uses the deferred native-modulo helper <see cref="ULongExtensions.MulMod64GpuCompatibleDeferred(ulong, ulong, ulong)"/>;
|
/// Uses the deferred native-modulo helper <see cref="ULongExtensions.MulMod64GpuDeferred(ulong, ulong, ulong)"/>;
|
||||||
/// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
|
/// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
|
||||||
/// </summary>
|
/// </summary>
|
||||||
/// <remarks>
|
/// <remarks>
|
||||||
@@ -186,7 +186,7 @@ public class MulMod64Benchmarks
|
|||||||
{
|
{
|
||||||
// TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks;
|
// 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.
|
// 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)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
|||||||
@@ -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);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
@@ -1,5 +1,6 @@
|
|||||||
using ILGPU;
|
using ILGPU;
|
||||||
using ILGPU.Runtime;
|
using ILGPU.Runtime;
|
||||||
|
using PerfectNumbers.Core;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core.Gpu;
|
namespace PerfectNumbers.Core.Gpu;
|
||||||
|
|
||||||
@@ -16,7 +17,7 @@ internal static class DivisorByDivisorKernels
|
|||||||
}
|
}
|
||||||
|
|
||||||
ulong exponent = exponents[index];
|
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<ulong> gaps, ArrayView<byte> deltas, byte modulus)
|
public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView<ulong> gaps, ArrayView<byte> deltas, byte modulus)
|
||||||
@@ -148,6 +149,6 @@ internal static class DivisorByDivisorKernels
|
|||||||
|
|
||||||
ulong exponent = exponents[index];
|
ulong exponent = exponents[index];
|
||||||
|
|
||||||
results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true);
|
results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponent, keepMontgomery: true);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@@ -1,5 +1,6 @@
|
|||||||
using ILGPU;
|
using ILGPU;
|
||||||
using ILGPU.Runtime;
|
using ILGPU.Runtime;
|
||||||
|
using PerfectNumbers.Core;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core.Gpu;
|
namespace PerfectNumbers.Core.Gpu;
|
||||||
|
|
||||||
@@ -12,6 +13,6 @@ internal static class Pow2MontgomeryKernels
|
|||||||
byte keepMontgomery,
|
byte keepMontgomery,
|
||||||
ArrayView1D<ulong, Stride1D.Dense> results)
|
ArrayView1D<ulong, Stride1D.Dense> results)
|
||||||
{
|
{
|
||||||
results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0);
|
results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponents[index], keepMontgomery != 0);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@@ -4,6 +4,7 @@ using System.Runtime.CompilerServices;
|
|||||||
using ILGPU;
|
using ILGPU;
|
||||||
using ILGPU.Algorithms.Random;
|
using ILGPU.Algorithms.Random;
|
||||||
using ILGPU.Runtime;
|
using ILGPU.Runtime;
|
||||||
|
using PerfectNumbers.Core;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core.Gpu;
|
namespace PerfectNumbers.Core.Gpu;
|
||||||
|
|
||||||
@@ -346,7 +347,7 @@ internal static partial class PrimeOrderGpuHeuristics
|
|||||||
}
|
}
|
||||||
|
|
||||||
ulong reduced = phi / factor;
|
ulong reduced = phi / factor;
|
||||||
if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL)
|
if (ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, reduced, keepMontgomery: false) == 1UL)
|
||||||
{
|
{
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
@@ -1130,7 +1131,7 @@ internal static partial class PrimeOrderGpuHeuristics
|
|||||||
|
|
||||||
private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
|
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)
|
private static ulong CalculateByDoublingKernel(ulong prime)
|
||||||
@@ -1155,7 +1156,7 @@ internal static partial class PrimeOrderGpuHeuristics
|
|||||||
internal static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
|
internal static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
|
||||||
{
|
{
|
||||||
ulong exponent = exponents[index];
|
ulong exponent = exponents[index];
|
||||||
remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false);
|
remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedGpu(divisor, exponent, keepMontgomery: false);
|
||||||
}
|
}
|
||||||
|
|
||||||
private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)
|
private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)
|
||||||
|
|||||||
@@ -3,6 +3,7 @@ using System.Collections.Concurrent;
|
|||||||
using ILGPU;
|
using ILGPU;
|
||||||
using ILGPU.Runtime;
|
using ILGPU.Runtime;
|
||||||
using System.Numerics;
|
using System.Numerics;
|
||||||
|
using PerfectNumbers.Core;
|
||||||
using System.Runtime.InteropServices;
|
using System.Runtime.InteropServices;
|
||||||
using System.Runtime.CompilerServices;
|
using System.Runtime.CompilerServices;
|
||||||
using ILGPU.Runtime.OpenCL;
|
using ILGPU.Runtime.OpenCL;
|
||||||
@@ -651,7 +652,7 @@ internal static partial class PrimeOrderGpuHeuristics
|
|||||||
return 0UL;
|
return 0UL;
|
||||||
}
|
}
|
||||||
|
|
||||||
return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisorData, keepMontgomery: false);
|
return ULongExtensions.Pow2MontgomeryModWindowedGpu(divisorData, exponent, keepMontgomery: false);
|
||||||
}
|
}
|
||||||
|
|
||||||
private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
|
private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
|
||||||
|
|||||||
@@ -1,4 +1,5 @@
|
|||||||
using System;
|
using System;
|
||||||
|
using System.Numerics;
|
||||||
using System.Runtime.CompilerServices;
|
using System.Runtime.CompilerServices;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core;
|
namespace PerfectNumbers.Core;
|
||||||
@@ -9,7 +10,7 @@ public static partial class ULongExtensions
|
|||||||
public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
|
public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
|
||||||
{
|
{
|
||||||
ulong modulus = divisor.Modulus;
|
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)
|
// if (exponent == 0UL)
|
||||||
// {
|
// {
|
||||||
// return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
|
// return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
|
||||||
@@ -22,8 +23,8 @@ public static partial class ULongExtensions
|
|||||||
// return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
|
// return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
|
||||||
// }
|
// }
|
||||||
|
|
||||||
int bitLength = GetPortableBitLength(exponent);
|
int bitLength = GetPortableBitLengthCpu(exponent);
|
||||||
int windowSize = GetWindowSize(bitLength);
|
int windowSize = GetWindowSizeCpu(bitLength);
|
||||||
int oddPowerCount = 1 << (windowSize - 1);
|
int oddPowerCount = 1 << (windowSize - 1);
|
||||||
|
|
||||||
ulong result = divisor.MontgomeryOne;
|
ulong result = divisor.MontgomeryOne;
|
||||||
@@ -219,4 +220,48 @@ public static partial class ULongExtensions
|
|||||||
order++;
|
order++;
|
||||||
return pow == 1UL;
|
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;
|
||||||
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
@@ -1,6 +1,7 @@
|
|||||||
using System.Collections.Concurrent;
|
using System.Collections.Concurrent;
|
||||||
using System.Runtime.CompilerServices;
|
using System.Runtime.CompilerServices;
|
||||||
using ILGPU;
|
using ILGPU;
|
||||||
|
using ILGPU.Algorithms;
|
||||||
using ILGPU.Runtime;
|
using ILGPU.Runtime;
|
||||||
using PerfectNumbers.Core.Gpu;
|
using PerfectNumbers.Core.Gpu;
|
||||||
|
|
||||||
@@ -17,7 +18,7 @@ public static partial class ULongExtensions
|
|||||||
}
|
}
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
[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,
|
// 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.
|
// 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)]
|
[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
|
// TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the
|
||||||
// 5-40× slowdown seen across real-world operand distributions.
|
// 5-40× slowdown seen across real-world operand distributions.
|
||||||
@@ -34,6 +35,147 @@ public static partial class ULongExtensions
|
|||||||
return state.MulModWithNativeModulo(b, modulus);
|
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<ulong> 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)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
|
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);
|
// return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
|
||||||
// }
|
// }
|
||||||
|
|
||||||
if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult))
|
return Pow2MontgomeryGpuExecutor.Execute(exponent, divisor, keepMontgomery);
|
||||||
{
|
|
||||||
return gpuResult;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery);
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
private static void InitializeStandardOddPowers(ulong modulus, Span<ulong> 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;
|
||||||
|
}
|
||||||
|
|
||||||
|
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)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
@@ -75,9 +302,9 @@ public static partial class ULongExtensions
|
|||||||
{
|
{
|
||||||
private static readonly ConcurrentDictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, byte, ArrayView1D<ulong, Stride1D.Dense>>> KernelCache = new();
|
private static readonly ConcurrentDictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, byte, ArrayView1D<ulong, Stride1D.Dense>>> 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);
|
GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true);
|
||||||
var execution = lease.EnterExecutionScope();
|
var execution = lease.EnterExecutionScope();
|
||||||
@@ -86,7 +313,7 @@ public static partial class ULongExtensions
|
|||||||
// Keep this commented out. It should never happen in production code.
|
// Keep this commented out. It should never happen in production code.
|
||||||
// if (accelerator.AcceleratorType == AcceleratorType.CPU)
|
// if (accelerator.AcceleratorType == AcceleratorType.CPU)
|
||||||
// {
|
// {
|
||||||
// return false;
|
// return 0UL;
|
||||||
// }
|
// }
|
||||||
|
|
||||||
var kernel = KernelCache.GetOrAdd(accelerator, static accel =>
|
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.
|
// Keep this commented. We don't want to catch any exceptions. All should crash the scanner.
|
||||||
// catch (AcceleratorException)
|
// catch (AcceleratorException)
|
||||||
// {
|
// {
|
||||||
// return false;
|
// return 0UL;
|
||||||
// }
|
// }
|
||||||
// catch (NotSupportedException)
|
// catch (NotSupportedException)
|
||||||
// {
|
// {
|
||||||
// return false;
|
// return 0UL;
|
||||||
// }
|
// }
|
||||||
// Intentionally avoid exception handling here; any accelerator failure should crash the scanner.
|
// Intentionally avoid exception handling here; any accelerator failure should crash the scanner.
|
||||||
return true;
|
return result;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@@ -498,96 +498,6 @@ public static partial class ULongExtensions
|
|||||||
return result;
|
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)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
|
private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
|
||||||
{
|
{
|
||||||
@@ -616,47 +526,6 @@ public static partial class ULongExtensions
|
|||||||
return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime);
|
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)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
|
public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
|
||||||
|
|||||||
Reference in new issue
Block a user