diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs
index 10a92a9..2907fe6 100644
--- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs
+++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs
@@ -42,7 +42,7 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi
// Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
// They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
- var exponentStepper = new ExponentRemainderStepper(cachedData);
+ var exponentStepper = new ExponentRemainderStepperCpu(cachedData);
var cycleStepper = new CycleRemainderStepper(divisorCycle);
diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
index 239684b..39e4dbe 100644
--- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
@@ -196,7 +196,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
// }
// else
{
- if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(
+ if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
candidate,
prime,
divisorData,
diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
similarity index 71%
rename from PerfectNumbers.Core/ExponentRemainderStepper.cs
rename to PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
index 9d6ad3f..df66c67 100644
--- a/PerfectNumbers.Core/ExponentRemainderStepper.cs
+++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
@@ -1,10 +1,10 @@
namespace PerfectNumbers.Core;
///
-/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates.
+/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path.
/// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order.
///
-internal struct ExponentRemainderStepper
+internal struct ExponentRemainderStepperCpu
{
private readonly MontgomeryDivisorData _divisor;
private readonly ulong _modulus;
@@ -13,7 +13,7 @@ internal struct ExponentRemainderStepper
private ulong _previousExponent;
private ulong _currentMontgomery;
- public ExponentRemainderStepper(in MontgomeryDivisorData divisor)
+ public ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor)
{
_divisor = divisor;
_modulus = divisor.Modulus;
@@ -81,28 +81,6 @@ internal struct ExponentRemainderStepper
return ReduceCurrent();
}
- ///
- /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence.
- ///
- /// First exponent evaluated on the GPU path. The same divisor metadata must be reused for future calls.
- /// The canonical residue produced by 2^exponent mod divisor.
- public ulong InitializeGpu(ulong exponent)
- {
- InitializeGpuState(exponent);
- return ReduceCurrent();
- }
-
- ///
- /// Initializes the GPU stepping state and reports whether the seed exponent produces unity.
- ///
- /// First exponent evaluated on the GPU path. Must match the divisor supplied during construction.
- /// when the seed exponent produces a Montgomery residue equal to one.
- public bool InitializeGpuIsUnity(ulong exponent)
- {
- InitializeGpuState(exponent);
- return _currentMontgomery == _montgomeryOne;
- }
-
///
/// Advances the CPU stepping state to the next exponent and reports whether it evaluates to unity.
/// Callers must seed the stepper through or and then
@@ -121,30 +99,6 @@ internal struct ExponentRemainderStepper
return _currentMontgomery == _montgomeryOne;
}
- ///
- /// Advances the GPU stepping state to the next exponent in a strictly ascending sequence that reuses the same divisor.
- /// Callers must invoke beforehand.
- ///
- /// Next exponent in the ascending sequence.
- /// The canonical residue produced by 2^exponent mod divisor.
- public ulong ComputeNextGpu(ulong exponent)
- {
- AdvanceGpu(exponent);
- return ReduceCurrent();
- }
-
- ///
- /// Advances the GPU stepping state to the next exponent and reports whether it evaluates to unity.
- /// Callers must seed the stepper through or .
- ///
- /// Next exponent in the ascending sequence.
- /// when the canonical residue equals one.
- public bool ComputeNextGpuIsUnity(ulong exponent)
- {
- AdvanceGpu(exponent);
- return _currentMontgomery == _montgomeryOne;
- }
-
///
/// Seeds the stepper with a precomputed Montgomery residue produced by the caller.
/// The provided exponent becomes the new baseline and all subsequent calls must continue with greater exponents.
@@ -169,14 +123,6 @@ internal struct ExponentRemainderStepper
/// Montgomery residue delta corresponding to the exponent gap.
/// Outputs whether the resulting residue equals Montgomery one.
/// when the advance succeeds.
- private void AdvanceGpu(ulong exponent)
- {
- ulong delta = exponent - _previousExponent;
- ulong multiplier = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, delta);
- _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
- _previousExponent = exponent;
- }
-
public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity)
{
// TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead
@@ -194,11 +140,5 @@ internal struct ExponentRemainderStepper
_previousExponent = exponent;
}
- private void InitializeGpuState(ulong exponent)
- {
- _currentMontgomery = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, exponent);
- _previousExponent = exponent;
- }
-
private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);
}
diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs
new file mode 100644
index 0000000..9d652f5
--- /dev/null
+++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs
@@ -0,0 +1,156 @@
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core.Gpu;
+
+///
+/// Tracks successive residues for a single divisor while kernels advance exponent candidates on the GPU path.
+/// Callers must seed the stepper once per divisor and then supply exponents in strictly ascending order.
+///
+internal struct ExponentRemainderStepperGpu
+{
+ private readonly MontgomeryDivisorData _divisor;
+ private readonly ulong _modulus;
+ private readonly ulong _nPrime;
+ private readonly ulong _montgomeryOne;
+ private ulong _previousExponent;
+ private ulong _currentMontgomery;
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public ExponentRemainderStepperGpu(in MontgomeryDivisorData divisor)
+ {
+ _divisor = divisor;
+ _modulus = divisor.Modulus;
+ _nPrime = divisor.NPrime;
+ _montgomeryOne = divisor.MontgomeryOne;
+ _previousExponent = 0UL;
+ _currentMontgomery = divisor.MontgomeryOne;
+ }
+
+ public ulong PreviousExponent
+ {
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ get => _previousExponent;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public bool MatchesDivisor(in MontgomeryDivisorData divisor)
+ {
+ return _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne;
+ }
+
+ ///
+ /// Clears the cached Montgomery residue so the stepper can be reused for another divisor.
+ /// Callers must reinitialize the stepper before consuming any residues after a reset.
+ ///
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public void Reset()
+ {
+ _previousExponent = 0UL;
+ _currentMontgomery = _montgomeryOne;
+ }
+
+ ///
+ /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence.
+ /// The returned value is the canonical residue corresponding to the provided exponent.
+ ///
+ /// First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata.
+ /// The canonical residue produced by 2^exponent mod divisor.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public ulong InitializeGpu(ulong exponent)
+ {
+ InitializeGpuState(exponent);
+ return ReduceCurrent();
+ }
+
+ ///
+ /// Initializes the GPU stepping state and reports whether the seed exponent already yields unity.
+ ///
+ /// First exponent evaluated for the current divisor. Must match the divisor passed to the constructor.
+ /// when the seed exponent produces a Montgomery residue equal to one.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public bool InitializeGpuIsUnity(ulong exponent)
+ {
+ InitializeGpuState(exponent);
+ return _currentMontgomery == _montgomeryOne;
+ }
+
+ ///
+ /// Advances the GPU stepping state to the next exponent in a strictly ascending sequence that shares the same divisor.
+ /// Callers must invoke beforehand.
+ ///
+ /// Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent.
+ /// The canonical residue produced by 2^exponent mod divisor.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public ulong ComputeNextGpu(ulong exponent)
+ {
+ AdvanceGpu(exponent);
+ return ReduceCurrent();
+ }
+
+ ///
+ /// Advances the GPU stepping state to the next exponent and reports whether it evaluates to unity.
+ /// Callers must seed the stepper through or and then
+ /// pass exponents in strictly ascending order.
+ ///
+ /// Next exponent in the ascending sequence.
+ /// when the canonical residue equals one.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public bool ComputeNextGpuIsUnity(ulong exponent)
+ {
+ AdvanceGpu(exponent);
+ return _currentMontgomery == _montgomeryOne;
+ }
+
+ ///
+ /// Seeds the stepper with a precomputed Montgomery residue produced by the caller.
+ /// The provided exponent becomes the new baseline and all subsequent calls must continue with greater exponents.
+ ///
+ /// Exponent associated with .
+ /// Residue in Montgomery form produced externally for the configured divisor.
+ /// Outputs whether the supplied residue equals Montgomery one.
+ /// when the residue is accepted and cached.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity)
+ {
+ _currentMontgomery = montgomeryResult;
+ _previousExponent = exponent;
+ isUnity = _currentMontgomery == _montgomeryOne;
+ return true;
+ }
+
+ ///
+ /// Advances the stepper using a caller-supplied Montgomery delta that represents the gap to the next exponent.
+ /// Callers must seed the stepper first, reuse the same divisor metadata, and continue feeding exponents in ascending order so the delta remains valid.
+ ///
+ /// Next exponent in the ascending sequence.
+ /// Montgomery residue delta corresponding to the exponent gap.
+ /// Outputs whether the resulting residue equals Montgomery one.
+ /// when the advance succeeds.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity)
+ {
+ _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime);
+ _previousExponent = exponent;
+ isUnity = _currentMontgomery == _montgomeryOne;
+ return true;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private void InitializeGpuState(ulong exponent)
+ {
+ _currentMontgomery = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, exponent);
+ _previousExponent = exponent;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private void AdvanceGpu(ulong exponent)
+ {
+ ulong delta = exponent - _previousExponent;
+ ulong multiplier = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, delta);
+ _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
+ _previousExponent = exponent;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);
+}
diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs
index 7d5e482..2f835b5 100644
--- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs
+++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs
@@ -41,37 +41,34 @@ internal static class DivisorByDivisorKernels
while (index >= 0)
{
- result = ULongExtensions.MulMod64Gpu(result, result, modulus);
+ ulong currentBit = (exponent >> index) & 1UL;
+ ulong squared = MultiplyMod(result, result, modulus);
+ bool processWindow = currentBit != 0UL;
- if (((exponent >> index) & 1UL) == 0UL)
+ result = processWindow ? result : squared;
+ index -= (int)(currentBit ^ 1UL);
+
+ int windowStartCandidate = index - windowSize + 1;
+ int negativeMask = windowStartCandidate >> 31;
+ windowStartCandidate &= ~negativeMask;
+
+ int windowStart = processWindow ? GetNextSetBitIndex(exponent, windowStartCandidate) : windowStartCandidate;
+ int windowLength = processWindow ? index - windowStart + 1 : 0;
+
+ for (int square = 0; square < windowLength; square++)
+ {
+ result = MultiplyMod(result, result, modulus);
+ }
+
+ if (!processWindow)
{
- 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 = 1; square < windowLength; square++)
- {
- result = ULongExtensions.MulMod64Gpu(result, result, modulus);
- }
-
ulong mask = (1UL << windowLength) - 1UL;
ulong windowValue = (exponent >> windowStart) & mask;
ulong multiplier = ComputeWindowedOddPower(windowValue, modulus);
- result = ULongExtensions.MulMod64Gpu(result, multiplier, modulus);
-
+ result = MultiplyMod(result, multiplier, modulus);
index = windowStart - 1;
}
@@ -82,25 +79,29 @@ internal static class DivisorByDivisorKernels
private static ulong ComputeWindowedOddPower(ulong windowValue, ulong modulus)
{
ulong baseValue = 2UL % modulus;
+ if (windowValue == 1UL)
+ {
+ return baseValue;
+ }
+
ulong result = baseValue;
ulong remaining = (windowValue - 1UL) >> 1;
- ulong squareBase = ULongExtensions.MulMod64Gpu(baseValue, baseValue, modulus);
+ ulong squareBase = MultiplyMod(baseValue, baseValue, modulus);
while (remaining != 0UL)
{
if ((remaining & 1UL) != 0UL)
{
- result = ULongExtensions.MulMod64Gpu(result, squareBase, modulus);
+ result = MultiplyMod(result, squareBase, modulus);
}
remaining >>= 1;
-
if (remaining == 0UL)
{
break;
}
- squareBase = ULongExtensions.MulMod64Gpu(squareBase, squareBase, modulus);
+ squareBase = MultiplyMod(squareBase, squareBase, modulus);
}
return result;
@@ -124,6 +125,65 @@ internal static class DivisorByDivisorKernels
window = bitLength <= 6 ? clamped : window;
return window;
}
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static int GetNextSetBitIndex(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)]
+ private static ulong MultiplyMod(ulong left, ulong right, ulong modulus)
+ {
+ if (left == 0UL || right == 0UL)
+ {
+ return 0UL;
+ }
+
+ ulong threshold = ulong.MaxValue / left;
+ if (right <= threshold)
+ {
+ ulong product = left * right;
+ return product % modulus;
+ }
+
+ return MultiplyModWide(left, right, modulus);
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong MultiplyModWide(ulong left, ulong right, ulong modulus)
+ {
+ ulong result = 0UL;
+ ulong addend = left;
+ ulong remaining = right;
+
+ while (remaining != 0UL)
+ {
+ if ((remaining & 1UL) != 0UL)
+ {
+ result += addend;
+ if (result >= modulus)
+ {
+ result -= modulus;
+ }
+ }
+
+ addend <<= 1;
+ if (addend >= modulus)
+ {
+ addend -= modulus;
+ }
+
+ remaining >>= 1;
+ }
+
+ return result;
+ }
+
public static void EvaluateDivisorWithStepperKernel(
Index1D _,
MontgomeryDivisorData divisor,
@@ -140,7 +200,7 @@ internal static class DivisorByDivisorKernels
// return;
// }
- var exponentStepper = new ExponentRemainderStepper(divisor);
+ var exponentStepper = new ExponentRemainderStepperGpu(divisor);
ulong firstPrime = exponents[0];
bool firstUnity = exponentStepper.InitializeGpuIsUnity(firstPrime);
// EvenPerfectBitScanner only schedules divisors with non-zero cycle lengths (equal to the tested prime),
diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
index 9c0f291..a0490be 100644
--- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
@@ -307,7 +307,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ResolveDivisorCycle(ulong divisor, ulong prime, in MontgomeryDivisorData divisorData)
{
- if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
+ if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentGpu(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
{
return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
}
@@ -318,7 +318,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool CandidatePassesHeuristics(ulong candidate, bool lastIsSeven)
{
- // Mirror the CPU heuristics so the GPU path rejects obvious composites before computing cycles.
+ // Mirror the CPU-side divisor screen so the GPU path rejects multiples of 3, 5, 7, and 11
+ // before resolving cycles, keeping cache lookups aligned with the persisted snapshot.
ulong remainder10 = candidate % 10UL;
bool accept10 = lastIsSeven
? (remainder10 == 3UL || remainder10 == 7UL || remainder10 == 9UL)
@@ -469,7 +470,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
if (cycle == 0UL)
{
- if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, firstPrime, cachedData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
+ if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentGpu(divisor, firstPrime, cachedData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
{
cycle = MersenneDivisorCycles.CalculateCycleLength(divisor, cachedData, skipPrimeOrderHeuristic: primeOrderFailed);
}
diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs
index c3fa2b1..7395e85 100644
--- a/PerfectNumbers.Core/MersenneDivisorCycles.cs
+++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs
@@ -81,7 +81,7 @@ public class MersenneDivisorCycles
return heuristicCycle == exponent;
}
- if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL)
+ if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL)
{
return cycleLength == exponent;
}
@@ -245,7 +245,88 @@ public class MersenneDivisorCycles
return order;
}
- public static bool TryCalculateCycleLengthForExponent(
+ public static bool TryCalculateCycleLengthForExponentCpu(
+ ulong divisor,
+ ulong exponent,
+ in MontgomeryDivisorData divisorData,
+ out ulong cycleLength,
+ out bool primeOrderFailed)
+ {
+ cycleLength = 0UL;
+ primeOrderFailed = false;
+
+ if ((divisor & (divisor - 1UL)) == 0UL)
+ {
+ cycleLength = 1UL;
+ return true;
+ }
+
+ if (divisor <= 3UL || exponent <= 1UL)
+ {
+ cycleLength = CalculateCycleLength(divisor, divisorData);
+ return true;
+ }
+
+ ulong computedOrder = PrimeOrderCalculator.Calculate(
+ divisor,
+ previousOrder: null,
+ divisorData,
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
+ if (computedOrder != 0UL)
+ {
+ cycleLength = computedOrder;
+ return true;
+ }
+
+ primeOrderFailed = true;
+
+ ulong phi = divisor - 1UL;
+ // A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without
+ // branching in the hot loop.
+ // if (phi == 0UL)
+ // {
+ // return false;
+ // }
+
+ int twoCount = BitOperations.TrailingZeroCount(phi);
+ ulong reducedPhi = phi >> twoCount;
+ if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
+ {
+ return false;
+ }
+
+ ulong k = reducedPhi / exponent;
+ // k == 0 would indicate reducedPhi < exponent, which the EvenPerfectBitScanner call flow prevents. Document the
+ // assumption instead of branching.
+ // if (k == 0UL)
+ // {
+ // return false;
+ // }
+
+ Dictionary factorCounts = ThreadStaticPools.RentFactorCountDictionary();
+ bool success = false;
+
+ if (twoCount > 0)
+ {
+ factorCounts[2UL] = twoCount;
+ }
+
+ bool exponentAccumulated = AccumulateFactors(exponent, factorCounts);
+ bool factorsAccumulated = exponentAccumulated && AccumulateFactors(k, factorCounts);
+
+ if (factorsAccumulated)
+ {
+ cycleLength = ReduceOrder(divisorData, phi, factorCounts);
+ success = true;
+ primeOrderFailed = false;
+ }
+
+ factorCounts.Clear();
+ ThreadStaticPools.ReturnFactorCountDictionary(factorCounts);
+ return success;
+ }
+ public static bool TryCalculateCycleLengthForExponentGpu(
ulong divisor,
ulong exponent,
in MontgomeryDivisorData divisorData,
@@ -482,7 +563,7 @@ public class MersenneDivisorCycles
ulong[]? heapCandidateArray = null;
bool[]? heapEvaluationArray = null;
- ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+ ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
ulong order = initialOrder;
for (int i = 0; i < count; i++)
@@ -513,7 +594,7 @@ public class MersenneDivisorCycles
ProcessReduceOrderPrime(heapCandidateArray.AsSpan(0, multiplicity), heapEvaluationArray.AsSpan(0, multiplicity), ref order, prime, multiplicity, ref stepper);
}
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
ThreadStaticPools.UlongPool.Return(primes);
if (heapCandidateArray is not null)
@@ -525,7 +606,7 @@ public class MersenneDivisorCycles
return order;
}
- private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepper stepper)
+ private static void ProcessReduceOrderPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper)
{
ulong working = order;
int actual = 0;
diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
index 5372881..ac93691 100644
--- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
+++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
@@ -163,12 +163,12 @@ internal static partial class PrimeOrderCalculator
Span candidates = buffer[..actual];
candidates.Sort();
- ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+ ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
int candidateCount = candidates.Length;
if (stepper.InitializeCpuIsUnity(candidates[0]))
{
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
return false;
}
@@ -176,12 +176,12 @@ internal static partial class PrimeOrderCalculator
{
if (stepper.ComputeNextIsUnity(candidates[i]))
{
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
return false;
}
}
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
return true;
}
@@ -219,7 +219,7 @@ internal static partial class PrimeOrderCalculator
buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value));
- ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+ ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
const int StackExponentCapacity = 16;
const int ExponentHardLimit = 256;
@@ -256,7 +256,7 @@ internal static partial class PrimeOrderCalculator
ProcessExponentLoweringPrime(heapCandidateArray.AsSpan(0, exponent), heapEvaluationArray.AsSpan(0, exponent), ref order, primeFactor, exponent, ref stepper);
}
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
if (heapCandidateArray is not null)
{
@@ -268,7 +268,7 @@ internal static partial class PrimeOrderCalculator
return order;
}
- private static void ProcessExponentLoweringPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepper stepper)
+ private static void ProcessExponentLoweringPrime(Span candidateBuffer, Span evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepperCpu stepper)
{
ulong working = order;
int actual = 0;
@@ -388,7 +388,7 @@ internal static partial class PrimeOrderCalculator
Span stackBuffer = stackalloc ulong[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
- ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+ ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool violates = false;
for (int i = 0; i < length; i++)
@@ -425,7 +425,7 @@ internal static partial class PrimeOrderCalculator
}
}
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
if (heapCandidateArray is not null)
{
@@ -435,7 +435,7 @@ internal static partial class PrimeOrderCalculator
return !violates;
}
- private static bool ValidateOrderForFactor(Span buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepper stepper)
+ private static bool ValidateOrderForFactor(Span buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepperCpu stepper)
{
ulong working = order;
int actual = 0;
@@ -547,7 +547,7 @@ internal static partial class PrimeOrderCalculator
int candidateCount = candidates.Count;
bool allowGpuBatch = true;
Span candidateSpan = CollectionsMarshal.AsSpan(candidates);
- ExponentRemainderStepper powStepper = ThreadStaticPools.RentExponentStepper(divisorData);
+ ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool powStepperInitialized = false;
// DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
@@ -626,7 +626,7 @@ internal static partial class PrimeOrderCalculator
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
- ThreadStaticPools.ReturnExponentStepper(powStepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
@@ -666,7 +666,7 @@ internal static partial class PrimeOrderCalculator
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
- ThreadStaticPools.ReturnExponentStepper(powStepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
@@ -681,7 +681,7 @@ internal static partial class PrimeOrderCalculator
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
- ThreadStaticPools.ReturnExponentStepper(powStepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
// DebugLog("No candidate confirmed");
return false;
}
@@ -892,7 +892,7 @@ internal static partial class PrimeOrderCalculator
Span stackBuffer = stackalloc ulong[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
- ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+ ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool violates = false;
for (int i = 0; i < length; i++)
@@ -929,7 +929,7 @@ internal static partial class PrimeOrderCalculator
}
}
- ThreadStaticPools.ReturnExponentStepper(stepper);
+ ThreadStaticPools.ReturnExponentStepperCpu(stepper);
if (heapCandidateArray is not null)
{
@@ -944,7 +944,7 @@ internal static partial class PrimeOrderCalculator
}
}
- private static bool CheckCandidateViolation(Span buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepper stepper)
+ private static bool CheckCandidateViolation(Span buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepperCpu stepper)
{
ulong working = candidate;
int actual = 0;
diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs
index 2da421e..477f8cf 100644
--- a/PerfectNumbers.Core/ThreadStaticPools.cs
+++ b/PerfectNumbers.Core/ThreadStaticPools.cs
@@ -329,31 +329,31 @@ namespace PerfectNumbers.Core
}
[ThreadStatic]
- private static ExponentRemainderStepper _exponentRemainderStepper;
+ private static ExponentRemainderStepperCpu _exponentRemainderStepperCpu;
[ThreadStatic]
- private static bool _hasExponentRemainderStepper;
+ private static bool _hasExponentRemainderStepperCpu;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- internal static ExponentRemainderStepper RentExponentStepper(in MontgomeryDivisorData divisorData)
+ internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData)
{
- if (_hasExponentRemainderStepper && _exponentRemainderStepper.MatchesDivisor(divisorData))
+ if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData))
{
- ExponentRemainderStepper stepper = _exponentRemainderStepper;
- _hasExponentRemainderStepper = false;
+ ExponentRemainderStepperCpu stepper = _exponentRemainderStepperCpu;
+ _hasExponentRemainderStepperCpu = false;
stepper.Reset();
return stepper;
}
- return new ExponentRemainderStepper(divisorData);
+ return new ExponentRemainderStepperCpu(divisorData);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- internal static void ReturnExponentStepper(ExponentRemainderStepper stepper)
+ internal static void ReturnExponentStepperCpu(ExponentRemainderStepperCpu stepper)
{
stepper.Reset();
- _exponentRemainderStepper = stepper;
- _hasExponentRemainderStepper = true;
+ _exponentRemainderStepperCpu = stepper;
+ _hasExponentRemainderStepperCpu = true;
}
}
}