Merge pull request #214 from mlemanczyk/codex/analyze-cpu-paths-in-evenperfectbitscanner

Restore cycle remainder TODO comments
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-10-20 16:07:08 +02:00
commit c726d6de6b
2 files changed
+376 -174

No files matched your search

+184 -17
View File
@@ -65,19 +65,16 @@ public static class UInt128Extensions
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow2MontgomeryModWindowed(this UInt128 exponent, UInt128 modulus) public static UInt128 Pow2MontgomeryModWindowed(this ulong exponent, UInt128 modulus)
{ {
// TODO: Prototype a UInt128-native Pow2MontgomeryModWindowed path that matches the GPU implementation once a faster modular multiplication helper is available, so we can revisit removing the GPU struct dependency without regressing benchmarks. // EvenPerfectBitScanner keeps CPU exponents above Pow2WindowFallbackThreshold once p >= 138,000,000,
if (modulus == UInt128.One) // so production never routes the tiny-exponent single-bit fallback through this helper. Leave the
{ // guard commented out for targeted microbenchmarks.
return UInt128.Zero; // This helper executes only on CPU fallbacks (PrimeOrderCalculator.Shared and PrimeOrderGpuHeuristics
} // call it after the GPU heuristics decline), but we still use GpuUInt128 because its mutable Montgomery
// operations beat the pure-UInt128 prototype by roughly 12% in Pow2MontgomeryModBenchmarks.
if (exponent == UInt128.Zero) // TODO: Swap in a CPU-native Montgomery ladder once it surpasses the shared struct so GPU and CPU tuning
{ // can diverge without compromising either path.
return UInt128.One % modulus;
}
GpuUInt128 modulusGpu = new(modulus); GpuUInt128 modulusGpu = new(modulus);
GpuUInt128 baseValue = new(2UL); GpuUInt128 baseValue = new(2UL);
if (baseValue.CompareTo(modulusGpu) >= 0) if (baseValue.CompareTo(modulusGpu) >= 0)
@@ -85,11 +82,11 @@ public static class UInt128Extensions
baseValue.Sub(modulusGpu); baseValue.Sub(modulusGpu);
} }
if (ShouldUseSingleBit(exponent)) // if (exponent <= Pow2WindowFallbackThreshold)
{ // {
GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponent, modulusGpu, baseValue); // GpuUInt128 tinyResult = Pow2MontgomeryModSingleBit((UInt128)exponent, modulusGpu, baseValue);
return (UInt128)singleBitResult; // return (UInt128)tinyResult;
} // }
GpuUInt128 exponentGpu = new(exponent); GpuUInt128 exponentGpu = new(exponent);
int bitLength = exponentGpu.GetBitLength(); int bitLength = exponentGpu.GetBitLength();
@@ -139,6 +136,176 @@ public static class UInt128Extensions
return (UInt128)result; return (UInt128)result;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow2MontgomeryModWindowed(this UInt128 exponent, UInt128 modulus)
{
// TODO: Prototype a UInt128-native Pow2MontgomeryModWindowed path that matches the GPU implementation once a faster
// modular multiplication helper is available, so we can replace the shared GpuUInt128 usage on CPU fallbacks and let
// the GPU kernels continue using their dedicated ladder without cross-contaminating optimizations.
// All GPU kernels already route through Pow2MontgomeryModWindowedGpu, and Pow2MontgomeryModBenchmarks confirmed this
// windowed ladder stays the fastest CPU option with the current Montgomery helpers.
if (modulus == UInt128.One)
{
return UInt128.Zero;
}
if (exponent == UInt128.Zero)
{
return UInt128.One % modulus;
}
// This helper still executes on CPU-only call sites, so we lean on the shared GpuUInt128 ladder until
// the TODO above lands and gives the CPU path its own Montgomery struct.
GpuUInt128 modulusGpu = new(modulus);
GpuUInt128 baseValue = new(2UL);
if (baseValue.CompareTo(modulusGpu) >= 0)
{
baseValue.Sub(modulusGpu);
}
// EvenPerfectBitScanner CPU exponents exceed Pow2WindowFallbackThreshold once p >= 138,000,000,
// so the single-bit fallback never triggers during production scans. Leave it commented out for targeted benchmarks.
// if (ShouldUseSingleBit(exponent))
// {
// GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponent, modulusGpu, baseValue);
// return (UInt128)singleBitResult;
// }
GpuUInt128 exponentGpu = new(exponent);
int bitLength = exponentGpu.GetBitLength();
int windowSize = GetWindowSize(bitLength);
int oddPowerCount = 1 << (windowSize - 1);
Span<GpuUInt128> oddPowers = stackalloc GpuUInt128[oddPowerCount];
InitializeOddPowers(baseValue, modulusGpu, oddPowers);
GpuUInt128 result = GpuUInt128.One;
int index = bitLength - 1;
while (index >= 0)
{
if (!IsBitSet(exponentGpu, index))
{
result.MulMod(result, modulusGpu);
index--;
continue;
}
int windowStart = index - windowSize + 1;
if (windowStart < 0)
{
windowStart = 0;
}
while (!IsBitSet(exponentGpu, windowStart))
{
windowStart++;
}
int windowBitCount = index - windowStart + 1;
for (int square = 0; square < windowBitCount; square++)
{
result.MulMod(result, modulusGpu);
}
ulong windowValue = ExtractWindowValue(exponentGpu, windowStart, windowBitCount);
int tableIndex = (int)((windowValue - 1UL) >> 1);
GpuUInt128 factor = oddPowers[tableIndex];
result.MulMod(factor, modulusGpu);
index = windowStart - 1;
}
return (UInt128)result;
}
/// <summary>
/// Uses the known cycle length to reduce the exponent and reuse cached remainders before
/// running the windowed pow2 ladder.
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow2MontgomeryModWindowedWithCycle(this UInt128 exponent, UInt128 modulus, ulong cycleLength, UInt128? precomputedRemainder = null)
{
// The EvenPerfectBitScanner CPU path receives cycle lengths from MersenneDivisorCycles, which never
// produces a zero-length order for admissible divisors. Keep the fallback documented for tests.
// if (cycleLength == 0UL)
// {
// return exponent.Pow2MontgomeryModWindowed(modulus);
// }
UInt128 remainder;
if (precomputedRemainder.HasValue)
{
// PowModWithCycle caches the reduced rotation count before reaching this helper, so production
// scans always reuse the precomputed remainder instead of recomputing the modulus here.
remainder = precomputedRemainder.Value;
}
else if (cycleLength != 0UL && exponent >= cycleLength)
{
if ((exponent >> 64) == UInt128.Zero)
{
remainder = (UInt128)((ulong)exponent % cycleLength);
}
else
{
remainder = exponent % (UInt128)cycleLength;
}
}
else
{
remainder = exponent;
}
// When the cycle divides the exponent, 2^exponent ≡ 1 (mod modulus) and the caller expects an
// immediate success signal.
if (remainder == UInt128.Zero)
{
return UInt128.One;
}
return remainder.Pow2MontgomeryModWindowed(modulus);
}
/// <summary>
/// Uses the known cycle length to reduce the exponent and reuse cached remainders before
/// running the windowed pow2 ladder.
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow2MontgomeryModWindowedWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength, UInt128? precomputedRemainder = null)
{
// Divisor-cycle lookups never hand the EvenPerfectBitScanner CPU flow a zero cycle length. Keep the
// fallback documented for targeted benchmarks and tests.
// if (cycleLength == UInt128.Zero)
// {
// return exponent.Pow2MontgomeryModWindowed(modulus);
// }
UInt128 remainder;
if (precomputedRemainder.HasValue)
{
// PowModWithCycle forwards the cached rotation count, so production scans always reuse the
// supplied remainder instead of recomputing exponent % cycleLength.
remainder = precomputedRemainder.Value;
}
else if (cycleLength != UInt128.Zero && exponent >= cycleLength)
{
remainder = exponent % cycleLength;
}
else
{
remainder = exponent;
}
// A zero remainder means the divisor order divides the exponent, so 2^exponent ≡ 1 (mod modulus)
// and the sieves rely on the early exit.
if (remainder == UInt128.Zero)
{
return UInt128.One;
}
return remainder.Pow2MontgomeryModWindowed(modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool ShouldUseSingleBit(UInt128 exponent) => (exponent >> 64) == UInt128.Zero && (ulong)exponent <= Pow2WindowFallbackThreshold; private static bool ShouldUseSingleBit(UInt128 exponent) => (exponent >> 64) == UInt128.Zero && (ulong)exponent <= Pow2WindowFallbackThreshold;
+192 -157
View File
@@ -10,6 +10,8 @@ public static partial class ULongExtensions
private const ulong Pow2WindowFallbackThreshold = 32UL; private const ulong Pow2WindowFallbackThreshold = 32UL;
private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1);
private static readonly UInt128 Pow2WindowedModulusThreshold = (UInt128)PerfectNumberConstants.MaxQForDivisorCycles;
private static readonly byte[] BitStatsZeroCount = new byte[256]; private static readonly byte[] BitStatsZeroCount = new byte[256];
private static readonly byte[] BitStatsPrefixZero = new byte[256]; private static readonly byte[] BitStatsPrefixZero = new byte[256];
private static readonly byte[] BitStatsSuffixZero = new byte[256]; private static readonly byte[] BitStatsSuffixZero = new byte[256];
@@ -593,218 +595,251 @@ public static partial class ULongExtensions
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 PowMod(this ulong exponent, UInt128 modulus) public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
{ {
UInt128 result = UInt128.One; // EvenPerfectBitScanner starts CPU scans at p >= 138,000,000 (Program.cs),
ulong exponentLoopIndex = 0UL; // so production never routes exponent == 0 here. Keep the guard commented out
// for targeted micro-benchmarks.
// if (exponent == 0UL)
// {
// return UInt128.One;
// }
// TODO: Port this scalar PowMod fallback to the ProcessEightBitWindows helper so CPU callers get the // MersenneResidueAutomaton seeds q = 2 * p * k + 1 (ModuloAutomata.cs),
// eight-bit window wins measured against the classic square-and-subtract implementation. // making modulus >= 2 * p + 1 > 276,000,000 and odd on the CPU path. The
// Return 1 because 2^0 = 1 // modulus == 1 guard is therefore unreachable during EvenPerfectBitScanner
if (exponent == 0UL) // runs.
return result; // if (modulus == UInt128.One)
// {
// return UInt128.Zero;
// }
// Any number mod 1 is 0 // Production exponents stay above 64 bits and modulus never drops below 4
if (modulus == UInt128.One) // for admissible divisors, so the classic doubling fallback does not
return UInt128.Zero; // trigger in scans.
// if (exponent < 64 || modulus < 4)
// {
// return Pow2ByDoubling(exponent, modulus);
// }
// For small exponents, do classic method // Every admissible divisor q is odd (2 * p * k + 1), preventing the
if (exponent < 64 || modulus < 4) // power-of-two mask path from activating in production.
// if ((modulus & (modulus - 1)) == 0)
// {
// UInt128 masked = UInt128.One << (int)(exponent & 127);
// return masked & (modulus - 1);
// }
if (!ShouldUseWindowedPow2(modulus))
{ {
for (; exponentLoopIndex < exponent; exponentLoopIndex++) return Pow2ByDoubling(exponent, modulus);
{
result <<= 1;
if (result >= modulus)
result -= modulus;
} }
return result; return exponent.Pow2MontgomeryModWindowed(modulus);
} }
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - 1)) == 0)
{
result = UInt128.One << (int)(exponent & 127);
return result & (modulus - 1);
}
// Reusing exponentLoopIndex to iterate again over the exponent range for the general-case accumulation.
exponentLoopIndex = 0UL;
// Reusing result after resetting it for the general modulus accumulation phase.
result = UInt128.One;
for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{
result <<= 1;
if (result >= modulus)
result -= modulus;
}
return result;
}
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using a known cycle length. /// Computes 2^exponent mod modulus using a known cycle length.
/// </summary> /// </summary>
public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength) public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength)
{ {
UInt128 result = UInt128.One; // EvenPerfectBitScanner keeps exponent >= 138,000,000 on CPU (Program.cs),
ulong exponentLoopIndex = 0UL; // so the zero-exponent guard stays disabled outside tests.
// if (exponent == 0UL)
// {
// return UInt128.One;
// }
// TODO: Wire this cycle-aware overload into the ProcessEightBitWindows helper so the reduced exponent path // The divisor generator produces odd q = 2 * p * k + 1 (ModuloAutomata.cs),
// inherits the faster windowed pow2 routine highlighted in the Pow2Montgomery benchmarks. // so modulus never reaches 1 in production.
// Return 1 because 2^0 = 1 // if (modulus == UInt128.One)
if (exponent == 0UL) // {
return result; // return UInt128.Zero;
// }
// Any number mod 1 is 0 // Every admissible divisor is odd, so the power-of-two shortcut cannot
if (modulus == UInt128.One) // trigger while scanning real candidates.
return UInt128.Zero; // if ((modulus & (modulus - 1)) == 0)
// {
// UInt128 masked = UInt128.One << (int)(exponent & 127);
// return masked & (modulus - 1);
// }
// For small exponents, do classic method bool hasCycle = cycleLength != 0UL;
if (exponent < 64 || modulus < 4)
{
for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{
result <<= 1;
if (result >= modulus)
result -= modulus;
}
return result;
}
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - 1)) == 0)
{
result = UInt128.One << (int)(exponent & 127);
return result & (modulus - 1);
}
// Reusing exponentLoopIndex to iterate over the rotation count for the cycle-aware path.
exponentLoopIndex = 0UL;
// Reusing result after resetting it for the rotation-accumulation pass.
result = UInt128.One;
// TODO: Replace this modulo with the cached cycle remainder produced by the divisor-cycle cache so PowModWithCycle avoids // TODO: Replace this modulo with the cached cycle remainder produced by the divisor-cycle cache so PowModWithCycle avoids
// repeated `%` work, matching the ProcessEightBitWindows wins captured in Pow2MontgomeryModCycleComputationBenchmarks. // repeated `%` work, matching the ProcessEightBitWindows wins captured in Pow2MontgomeryModCycleComputationBenchmarks.
ulong rotationCount = exponent % cycleLength; ulong rotationCount = hasCycle && exponent >= cycleLength ? exponent % cycleLength : exponent;
for (; exponentLoopIndex < rotationCount; exponentLoopIndex++)
// Keep this fast path: valid divisors satisfy cycle | exponent, producing
// a zero remainder that means 2^exponent ≡ 1 (mod q).
if (hasCycle && rotationCount == 0UL)
{ {
result <<= 1; return UInt128.One;
if (result >= modulus)
result -= modulus;
} }
return result; // FastDiv64 can return tiny divisors (see MersenneNumberIncrementalCpuTester),
// so the small-rotation fallback still fires on the CPU path.
if (modulus < 4 || rotationCount < 64)
{
return Pow2ByDoubling(rotationCount, modulus);
} }
if (ShouldUseWindowedPow2(modulus))
{
UInt128 exponent128 = exponent;
UInt128 rotationCount128 = rotationCount;
return exponent128.Pow2MontgomeryModWindowedWithCycle(modulus, cycleLength, rotationCount128);
}
return Pow2ByDoubling(rotationCount, modulus);
}
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using a known cycle length. /// Computes 2^exponent mod modulus using a known cycle length.
/// </summary> /// </summary>
public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength) public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength)
{ {
UInt128 result = UInt128.One; // CPU scans only feed exponents >= 138,000,000 (Program.cs), so keep the
ulong exponentLoopIndex = 0UL; // zero-exponent shortcut commented out for tests.
// if (exponent == 0UL)
// {
// return UInt128.One;
// }
// TODO: Replace this UInt128-cycle overload with the ProcessEightBitWindows helper so large-exponent CPU scans // Modulus originates from q = 2 * p * k + 1 (ModuloAutomata.cs), therefore
// reuse the faster windowed pow2 ladder instead of the manual rotation loop measured to lag behind in benchmarks. // modulus == 1 never occurs outside dedicated benchmarks.
// Return 1 because 2^0 = 1 // if (modulus == UInt128.One)
if (exponent == UInt128.Zero) // {
return result; // return UInt128.Zero;
// }
// Any number mod 1 is 0 // Admissible q are odd, rendering the power-of-two shortcut unreachable.
if (modulus == UInt128.One) // if ((modulus & (modulus - UInt128.One)) == UInt128.Zero)
return UInt128.Zero; // {
// UInt128 masked = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
// return masked & (modulus - UInt128.One);
// }
// For small exponents, do classic method bool hasCycle = cycleLength != UInt128.Zero;
if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four) UInt128 exponent128 = exponent;
{
for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{
result <<= 1;
if (result >= modulus)
result -= modulus;
}
return result;
}
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - UInt128.One)) == UInt128.Zero)
{
result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
return result & (modulus - 1);
}
// Reusing result after resetting it for the rotation-driven accumulation phase.
result = UInt128.One;
// TODO: Swap this modulo with the upcoming UInt128 cycle remainder helper so large-exponent scans reuse cached // TODO: Swap this modulo with the upcoming UInt128 cycle remainder helper so large-exponent scans reuse cached
// reductions instead of recomputing `%` for every lookup, as demonstrated in Pow2MontgomeryModCycleComputationBenchmarks. // reductions instead of recomputing `%` for every lookup, as demonstrated in Pow2MontgomeryModCycleComputationBenchmarks.
UInt128 rotationCount = exponent % cycleLength; UInt128 rotationCount = hasCycle && exponent128 >= cycleLength ? exponent128 % cycleLength : exponent128;
UInt128 rotationIndex = UInt128.Zero;
while (rotationIndex < rotationCount) if (hasCycle && rotationCount == UInt128.Zero)
{ {
result <<= 1; return UInt128.One;
if (result >= modulus)
result -= modulus;
rotationIndex += UInt128.One;
} }
return result; if (modulus < UInt128Numbers.Four || rotationCount < UInt128Numbers.SixtyFour)
{
return Pow2ByDoubling(rotationCount, modulus);
} }
if (ShouldUseWindowedPow2(modulus))
{
return exponent128.Pow2MontgomeryModWindowedWithCycle(modulus, cycleLength, rotationCount);
}
return Pow2ByDoubling(rotationCount, modulus);
}
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using a known cycle length. /// Computes 2^exponent mod modulus using a known cycle length.
/// </summary> /// </summary>
public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength) public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength)
{ {
UInt128 one = UInt128.One, UInt128 one = UInt128.One;
result = one, UInt128 zero = UInt128.Zero;
zero = UInt128.Zero;
ulong exponentLoopIndex = 0UL;
// TODO: Migrate this UInt128 exponent overload to ProcessEightBitWindows so the large-cycle reductions drop the // EvenPerfectBitScanner drives this path with exponents derived from p, phi,
// slow manual loop that underperforms the windowed pow2 helper in the Pow2 benchmark suite. // or divisor quotients, none of which are zero on the production scan.
// Return 1 because 2^0 = 1 // if (exponent == zero)
if (exponent == zero) // {
return result; // return one;
// }
// Any number mod 1 is 0 // Modulus again comes from q = 2 * p * k + 1, so modulus == 1 is unreachable.
if (modulus == one) // if (modulus == one)
return zero; // {
// return zero;
// }
// For small exponents, do classic method // Candidate divisors stay odd, therefore the power-of-two shortcut stays off.
if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four) // if ((modulus & (modulus - one)) == zero)
{ // {
for (; exponentLoopIndex < exponent; exponentLoopIndex++) // UInt128 masked = one << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
{ // return masked & (modulus - one);
result <<= 1; // }
if (result >= modulus)
result -= modulus;
}
return result; bool hasCycle = cycleLength != zero;
}
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - one)) == zero)
{
result = one << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
return result & (modulus - 1);
}
// Reusing result after resetting it for the rotation-driven accumulation phase.
result = one;
// TODO: Swap this modulo with the shared UInt128 cycle remainder helper once available so CRT powmods reuse cached // TODO: Swap this modulo with the shared UInt128 cycle remainder helper once available so CRT powmods reuse cached
// reductions in the windowed ladder, avoiding the `%` cost highlighted in Pow2MontgomeryModCycleComputationBenchmarks. // reductions in the windowed ladder, avoiding the `%` cost highlighted in Pow2MontgomeryModCycleComputationBenchmarks.
UInt128 rotationCount = exponent % cycleLength; UInt128 rotationCount = hasCycle && exponent >= cycleLength ? exponent % cycleLength : exponent;
// We're reusing "zero" as rotation index for just a little better performance if (hasCycle && rotationCount == zero)
while (zero < rotationCount) {
return one;
}
if (modulus < UInt128Numbers.Four || rotationCount < UInt128Numbers.SixtyFour)
{
return Pow2ByDoubling(rotationCount, modulus);
}
if (ShouldUseWindowedPow2(modulus))
{
return exponent.Pow2MontgomeryModWindowedWithCycle(modulus, cycleLength, rotationCount);
}
return Pow2ByDoubling(rotationCount, modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool ShouldUseWindowedPow2(UInt128 modulus) => modulus > Pow2WindowedModulusThreshold;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static UInt128 Pow2ByDoubling(ulong exponent, UInt128 modulus)
{
UInt128 result = UInt128.One;
for (ulong i = 0; i < exponent; i++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) while (result >= modulus)
{
result -= modulus; result -= modulus;
}
}
zero += one; return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static UInt128 Pow2ByDoubling(UInt128 exponent, UInt128 modulus)
{
UInt128 result = UInt128.One;
UInt128 index = UInt128.Zero;
while (index < exponent)
{
result <<= 1;
while (result >= modulus)
{
result -= modulus;
}
index += UInt128.One;
} }
return result; return result;