BigInteger-friendly steppers added. MersenneResidueStepper converted to struct. Manual residue tracking replaced with steppers. Modulo operations replaced with steppers.
This commit is contained in:
1 parent
fbf074ade0
commit
6887111128
9 files changed
+208
-391
No files matched your search
@@ -0,0 +1,38 @@
|
|||||||
|
using System.Numerics;
|
||||||
|
using System.Runtime.CompilerServices;
|
||||||
|
|
||||||
|
namespace PerfectNumbers.Core;
|
||||||
|
|
||||||
|
/// <summary>
|
||||||
|
/// Maintains rolling remainders of an increasing sequence modulo a BigInteger cycle length.
|
||||||
|
/// Mirrors <see cref="CycleRemainderStepper"/> for values that exceed 64 bits.
|
||||||
|
/// </summary>
|
||||||
|
[method: MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public struct BigIntegerCycleRemainderStepper(BigInteger cycleLength)
|
||||||
|
{
|
||||||
|
private readonly BigInteger _cycleLength = cycleLength;
|
||||||
|
private BigInteger _previousValue = BigInteger.Zero;
|
||||||
|
private BigInteger _currentRemainder = BigInteger.Zero;
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public BigInteger Initialize(BigInteger value)
|
||||||
|
{
|
||||||
|
_previousValue = value;
|
||||||
|
_currentRemainder = value % _cycleLength;
|
||||||
|
return _currentRemainder;
|
||||||
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public BigInteger ComputeNext(BigInteger value)
|
||||||
|
{
|
||||||
|
BigInteger delta = value - _previousValue;
|
||||||
|
_previousValue = value;
|
||||||
|
_currentRemainder += delta;
|
||||||
|
if (_currentRemainder >= _cycleLength)
|
||||||
|
{
|
||||||
|
_currentRemainder %= _cycleLength;
|
||||||
|
}
|
||||||
|
|
||||||
|
return _currentRemainder;
|
||||||
|
}
|
||||||
|
}
|
||||||
@@ -0,0 +1,60 @@
|
|||||||
|
using System.Numerics;
|
||||||
|
using System.Runtime.CompilerServices;
|
||||||
|
|
||||||
|
namespace PerfectNumbers.Core;
|
||||||
|
|
||||||
|
/// <summary>
|
||||||
|
/// Tracks successive residues of <c>baseValue^exponent mod modulus</c> for ascending exponents using <see cref="BigInteger"/> inputs.
|
||||||
|
/// This avoids recomputing full modular exponentiation when stepping forward by a known delta.
|
||||||
|
/// </summary>
|
||||||
|
public struct BigIntegerExponentStepper(BigInteger modulus, BigInteger baseValue)
|
||||||
|
{
|
||||||
|
private readonly BigInteger _modulus = modulus;
|
||||||
|
private readonly BigInteger _baseValue = baseValue % modulus;
|
||||||
|
private BigInteger PreviousExponent = BigInteger.Zero;
|
||||||
|
private BigInteger CurrentResidue = BigInteger.One % modulus;
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public BigInteger Initialize(BigInteger exponent)
|
||||||
|
{
|
||||||
|
CurrentResidue = BigInteger.ModPow(_baseValue, exponent, _modulus);
|
||||||
|
PreviousExponent = exponent;
|
||||||
|
return CurrentResidue;
|
||||||
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public bool InitializeIsUnity(BigInteger exponent)
|
||||||
|
{
|
||||||
|
Initialize(exponent);
|
||||||
|
return CurrentResidue.IsOne;
|
||||||
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public BigInteger ComputeNext(BigInteger exponent)
|
||||||
|
{
|
||||||
|
if (exponent == PreviousExponent)
|
||||||
|
{
|
||||||
|
return CurrentResidue;
|
||||||
|
}
|
||||||
|
|
||||||
|
BigInteger delta = exponent - PreviousExponent;
|
||||||
|
BigInteger deltaResidue = BigInteger.ModPow(_baseValue, delta, _modulus);
|
||||||
|
CurrentResidue = (CurrentResidue * deltaResidue) % _modulus;
|
||||||
|
PreviousExponent = exponent;
|
||||||
|
return CurrentResidue;
|
||||||
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public bool ComputeNextIsUnity(BigInteger exponent)
|
||||||
|
{
|
||||||
|
ComputeNext(exponent);
|
||||||
|
return CurrentResidue.IsOne;
|
||||||
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public void Reset()
|
||||||
|
{
|
||||||
|
PreviousExponent = BigInteger.Zero;
|
||||||
|
CurrentResidue = BigInteger.One % _modulus;
|
||||||
|
}
|
||||||
|
}
|
||||||
@@ -2,6 +2,7 @@ using System.Numerics;
|
|||||||
using System.Runtime.CompilerServices;
|
using System.Runtime.CompilerServices;
|
||||||
using System.Globalization;
|
using System.Globalization;
|
||||||
using System.IO;
|
using System.IO;
|
||||||
|
using PerfectNumbers.Core;
|
||||||
using PerfectNumbers.Core.Gpu;
|
using PerfectNumbers.Core.Gpu;
|
||||||
using PerfectNumbers.Core.Gpu.Accelerators;
|
using PerfectNumbers.Core.Gpu.Accelerators;
|
||||||
|
|
||||||
@@ -294,40 +295,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
ulong limitHigh = limit.High;
|
ulong limitHigh = limit.High;
|
||||||
ulong divisorLow = divisor.Low;
|
ulong divisorLow = divisor.Low;
|
||||||
|
|
||||||
// Keep the divisibility filters aligned with the divisor-cycle generator so the
|
|
||||||
// CPU path never requests cycles that were skipped during cache creation.
|
|
||||||
LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
|
|
||||||
ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);
|
|
||||||
|
|
||||||
byte step10;
|
|
||||||
byte step8;
|
|
||||||
byte step3;
|
|
||||||
byte step7;
|
|
||||||
byte step11;
|
|
||||||
byte step13;
|
|
||||||
byte step17;
|
|
||||||
byte step19;
|
|
||||||
|
|
||||||
byte remainder10;
|
|
||||||
byte remainder8;
|
|
||||||
byte remainder3;
|
|
||||||
byte remainder7;
|
|
||||||
byte remainder11;
|
|
||||||
byte remainder13;
|
|
||||||
byte remainder17;
|
|
||||||
byte remainder19;
|
|
||||||
|
|
||||||
if (stepHigh == 0UL && limitHigh == 0UL)
|
if (stepHigh == 0UL && limitHigh == 0UL)
|
||||||
{
|
{
|
||||||
step10 = (byte)(stepLow % 10UL);
|
|
||||||
step8 = (byte)(stepLow & 7UL);
|
|
||||||
step3 = (byte)(stepLow % 3UL);
|
|
||||||
step7 = (byte)(stepLow % 7UL);
|
|
||||||
step11 = (byte)(stepLow % 11UL);
|
|
||||||
step13 = (byte)(stepLow % 13UL);
|
|
||||||
step17 = (byte)(stepLow % 17UL);
|
|
||||||
step19 = (byte)(stepLow % 19UL);
|
|
||||||
|
|
||||||
ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL;
|
ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL;
|
||||||
if (maxK == 0UL)
|
if (maxK == 0UL)
|
||||||
{
|
{
|
||||||
@@ -350,15 +319,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
allowedMax,
|
allowedMax,
|
||||||
startK,
|
startK,
|
||||||
maxK,
|
maxK,
|
||||||
decimalMask,
|
|
||||||
step10,
|
|
||||||
step8,
|
|
||||||
step3,
|
|
||||||
step7,
|
|
||||||
step11,
|
|
||||||
step13,
|
|
||||||
step17,
|
|
||||||
step19,
|
|
||||||
ref currentK,
|
ref currentK,
|
||||||
out processedTop,
|
out processedTop,
|
||||||
out foundDivisor);
|
out foundDivisor);
|
||||||
@@ -381,15 +341,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
allowedMax,
|
allowedMax,
|
||||||
1UL,
|
1UL,
|
||||||
lowerEnd,
|
lowerEnd,
|
||||||
decimalMask,
|
|
||||||
step10,
|
|
||||||
step8,
|
|
||||||
step3,
|
|
||||||
step7,
|
|
||||||
step11,
|
|
||||||
step13,
|
|
||||||
step17,
|
|
||||||
step19,
|
|
||||||
ref currentK,
|
ref currentK,
|
||||||
out processedBottom,
|
out processedBottom,
|
||||||
out foundDivisor);
|
out foundDivisor);
|
||||||
@@ -400,29 +351,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
}
|
}
|
||||||
|
|
||||||
ulong divisorHigh = divisor.High;
|
ulong divisorHigh = divisor.High;
|
||||||
step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
|
var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor);
|
||||||
step8 = (byte)(stepLow % 8UL);
|
|
||||||
step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL);
|
|
||||||
step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL);
|
|
||||||
step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL);
|
|
||||||
step13 = (byte)((((stepHigh % 13UL) * 3UL) + (stepLow % 13UL)) % 13UL);
|
|
||||||
step17 = (byte)(((stepHigh % 17UL) + (stepLow % 17UL)) % 17UL);
|
|
||||||
step19 = (byte)((((stepHigh % 19UL) * 17UL) + (stepLow % 19UL)) % 19UL);
|
|
||||||
|
|
||||||
remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL);
|
|
||||||
remainder8 = (byte)(divisorLow % 8UL);
|
|
||||||
remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL);
|
|
||||||
remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
|
|
||||||
remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL);
|
|
||||||
remainder13 = (byte)((((divisorHigh % 13UL) * 3UL) + (divisorLow % 13UL)) % 13UL);
|
|
||||||
remainder17 = (byte)(((divisorHigh % 17UL) + (divisorLow % 17UL)) % 17UL);
|
|
||||||
remainder19 = (byte)((((divisorHigh % 19UL) * 17UL) + (divisorLow % 19UL)) % 19UL);
|
|
||||||
|
|
||||||
var divisorPool = MontgomeryDivisorDataPool.Shared;
|
var divisorPool = MontgomeryDivisorDataPool.Shared;
|
||||||
while (divisor.CompareTo(limit) <= 0)
|
while (divisor.CompareTo(limit) <= 0)
|
||||||
{
|
{
|
||||||
bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0;
|
if (residueStepper.IsAdmissible())
|
||||||
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
|
|
||||||
{
|
{
|
||||||
ulong candidate = divisor.Low;
|
ulong candidate = divisor.Low;
|
||||||
MontgomeryDivisorData divisorData = divisorPool.FromModulus(candidate);
|
MontgomeryDivisorData divisorData = divisorPool.FromModulus(candidate);
|
||||||
@@ -468,14 +402,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
|
|
||||||
divisor.Add(step);
|
divisor.Add(step);
|
||||||
currentK++;
|
currentK++;
|
||||||
remainder10 = AddMod10(remainder10, step10);
|
residueStepper.Advance();
|
||||||
remainder8 = AddMod8(remainder8, step8);
|
|
||||||
remainder3 = AddMod3(remainder3, step3);
|
|
||||||
remainder7 = AddMod7(remainder7, step7);
|
|
||||||
remainder11 = AddMod11(remainder11, step11);
|
|
||||||
remainder13 = AddMod13(remainder13, step13);
|
|
||||||
remainder17 = AddMod17(remainder17, step17);
|
|
||||||
remainder19 = AddMod19(remainder19, step19);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
processedAll = true;
|
processedAll = true;
|
||||||
@@ -510,27 +437,29 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
|
LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
|
||||||
ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);
|
ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);
|
||||||
|
|
||||||
byte step10 = (byte)(step % 10);
|
var rem10 = new BigIntegerCycleRemainderStepper(10);
|
||||||
byte step8 = (byte)(step % 8);
|
var rem8 = new BigIntegerCycleRemainderStepper(8);
|
||||||
byte step3 = (byte)(step % 3);
|
var rem5 = new BigIntegerCycleRemainderStepper(5);
|
||||||
byte step7 = (byte)(step % 7);
|
var rem3 = new BigIntegerCycleRemainderStepper(3);
|
||||||
byte step11 = (byte)(step % 11);
|
var rem7 = new BigIntegerCycleRemainderStepper(7);
|
||||||
byte step13 = (byte)(step % 13);
|
var rem11 = new BigIntegerCycleRemainderStepper(11);
|
||||||
byte step17 = (byte)(step % 17);
|
var rem13 = new BigIntegerCycleRemainderStepper(13);
|
||||||
byte step19 = (byte)(step % 19);
|
var rem17 = new BigIntegerCycleRemainderStepper(17);
|
||||||
|
var rem19 = new BigIntegerCycleRemainderStepper(19);
|
||||||
|
|
||||||
byte remainder10 = (byte)(divisor % 10);
|
byte remainder10 = (byte)rem10.Initialize(divisor);
|
||||||
byte remainder8 = (byte)(divisor % 8);
|
byte remainder8 = (byte)rem8.Initialize(divisor);
|
||||||
byte remainder3 = (byte)(divisor % 3);
|
byte remainder5 = (byte)rem5.Initialize(divisor);
|
||||||
byte remainder7 = (byte)(divisor % 7);
|
byte remainder3 = (byte)rem3.Initialize(divisor);
|
||||||
byte remainder11 = (byte)(divisor % 11);
|
byte remainder7 = (byte)rem7.Initialize(divisor);
|
||||||
byte remainder13 = (byte)(divisor % 13);
|
byte remainder11 = (byte)rem11.Initialize(divisor);
|
||||||
byte remainder17 = (byte)(divisor % 17);
|
byte remainder13 = (byte)rem13.Initialize(divisor);
|
||||||
byte remainder19 = (byte)(divisor % 19);
|
byte remainder17 = (byte)rem17.Initialize(divisor);
|
||||||
|
byte remainder19 = (byte)rem19.Initialize(divisor);
|
||||||
|
|
||||||
while (divisor <= allowedMax)
|
while (divisor <= allowedMax)
|
||||||
{
|
{
|
||||||
bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0;
|
bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0;
|
||||||
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
|
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
|
||||||
{
|
{
|
||||||
if (IsProbablePrimeBigInteger(divisor))
|
if (IsProbablePrimeBigInteger(divisor))
|
||||||
@@ -552,14 +481,15 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
|
|
||||||
currentK += BigInteger.One;
|
currentK += BigInteger.One;
|
||||||
divisor += step;
|
divisor += step;
|
||||||
remainder10 = AddMod10(remainder10, step10);
|
remainder10 = (byte)rem10.ComputeNext(divisor);
|
||||||
remainder8 = AddMod8(remainder8, step8);
|
remainder8 = (byte)rem8.ComputeNext(divisor);
|
||||||
remainder3 = AddMod3(remainder3, step3);
|
remainder5 = (byte)rem5.ComputeNext(divisor);
|
||||||
remainder7 = AddMod7(remainder7, step7);
|
remainder3 = (byte)rem3.ComputeNext(divisor);
|
||||||
remainder11 = AddMod11(remainder11, step11);
|
remainder7 = (byte)rem7.ComputeNext(divisor);
|
||||||
remainder13 = AddMod13(remainder13, step13);
|
remainder11 = (byte)rem11.ComputeNext(divisor);
|
||||||
remainder17 = AddMod17(remainder17, step17);
|
remainder13 = (byte)rem13.ComputeNext(divisor);
|
||||||
remainder19 = AddMod19(remainder19, step19);
|
remainder17 = (byte)rem17.ComputeNext(divisor);
|
||||||
|
remainder19 = (byte)rem19.ComputeNext(divisor);
|
||||||
|
|
||||||
if (divisor > allowedMax)
|
if (divisor > allowedMax)
|
||||||
{
|
{
|
||||||
@@ -579,15 +509,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
ulong allowedMax,
|
ulong allowedMax,
|
||||||
ulong startK,
|
ulong startK,
|
||||||
ulong endK,
|
ulong endK,
|
||||||
ushort decimalMask,
|
|
||||||
byte step10,
|
|
||||||
byte step8,
|
|
||||||
byte step3,
|
|
||||||
byte step7,
|
|
||||||
byte step11,
|
|
||||||
byte step13,
|
|
||||||
byte step17,
|
|
||||||
byte step19,
|
|
||||||
ref ulong currentK,
|
ref ulong currentK,
|
||||||
out bool processedAll,
|
out bool processedAll,
|
||||||
out ulong foundDivisor)
|
out ulong foundDivisor)
|
||||||
@@ -619,14 +540,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
ulong startDivisor = (ulong)startDivisor128;
|
ulong startDivisor = (ulong)startDivisor128;
|
||||||
ulong rangeLimit = (ulong)rangeLimit128;
|
ulong rangeLimit = (ulong)rangeLimit128;
|
||||||
|
|
||||||
byte remainder10 = (byte)(startDivisor % 10UL);
|
var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)step128, (GpuUInt128)startDivisor128);
|
||||||
byte remainder8 = (byte)(startDivisor & 7UL);
|
|
||||||
byte remainder3 = (byte)(startDivisor % 3UL);
|
|
||||||
byte remainder7 = (byte)(startDivisor % 7UL);
|
|
||||||
byte remainder11 = (byte)(startDivisor % 11UL);
|
|
||||||
byte remainder13 = (byte)(startDivisor % 13UL);
|
|
||||||
byte remainder17 = (byte)(startDivisor % 17UL);
|
|
||||||
byte remainder19 = (byte)(startDivisor % 19UL);
|
|
||||||
|
|
||||||
return CheckDivisors64(
|
return CheckDivisors64(
|
||||||
gpu,
|
gpu,
|
||||||
@@ -634,23 +548,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
|
|||||||
step,
|
step,
|
||||||
rangeLimit,
|
rangeLimit,
|
||||||
startDivisor,
|
startDivisor,
|
||||||
decimalMask,
|
ref residueStepper,
|
||||||
step10,
|
|
||||||
step8,
|
|
||||||
step3,
|
|
||||||
step7,
|
|
||||||
step11,
|
|
||||||
step13,
|
|
||||||
step17,
|
|
||||||
step19,
|
|
||||||
remainder10,
|
|
||||||
remainder8,
|
|
||||||
remainder3,
|
|
||||||
remainder7,
|
|
||||||
remainder11,
|
|
||||||
remainder13,
|
|
||||||
remainder17,
|
|
||||||
remainder19,
|
|
||||||
ref currentK,
|
ref currentK,
|
||||||
out processedAll,
|
out processedAll,
|
||||||
out foundDivisor);
|
out foundDivisor);
|
||||||
@@ -662,23 +560,7 @@ private bool CheckDivisors64(
|
|||||||
ulong step,
|
ulong step,
|
||||||
ulong limit,
|
ulong limit,
|
||||||
ulong divisor,
|
ulong divisor,
|
||||||
ushort decimalMask,
|
ref MersenneDivisorResidueStepper residueStepper,
|
||||||
byte step10,
|
|
||||||
byte step8,
|
|
||||||
byte step3,
|
|
||||||
byte step7,
|
|
||||||
byte step11,
|
|
||||||
byte step13,
|
|
||||||
byte step17,
|
|
||||||
byte step19,
|
|
||||||
byte remainder10,
|
|
||||||
byte remainder8,
|
|
||||||
byte remainder3,
|
|
||||||
byte remainder7,
|
|
||||||
byte remainder11,
|
|
||||||
byte remainder13,
|
|
||||||
byte remainder17,
|
|
||||||
byte remainder19,
|
|
||||||
ref ulong currentK,
|
ref ulong currentK,
|
||||||
out bool processedAll,
|
out bool processedAll,
|
||||||
out ulong foundDivisor)
|
out ulong foundDivisor)
|
||||||
@@ -687,9 +569,7 @@ private bool CheckDivisors64(
|
|||||||
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
|
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
|
||||||
if (!canAdvance)
|
if (!canAdvance)
|
||||||
{
|
{
|
||||||
if (divisor <= limit)
|
if (divisor <= limit && residueStepper.IsAdmissible())
|
||||||
{
|
|
||||||
if (remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0 && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
|
|
||||||
{
|
{
|
||||||
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
|
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
|
||||||
ulong divisorCycle;
|
ulong divisorCycle;
|
||||||
@@ -726,7 +606,6 @@ private bool CheckDivisors64(
|
|||||||
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
|
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
|
||||||
|
|
||||||
processedAll = true;
|
processedAll = true;
|
||||||
foundDivisor = 0UL;
|
foundDivisor = 0UL;
|
||||||
@@ -743,8 +622,7 @@ private bool CheckDivisors64(
|
|||||||
ulong remainingIterations = ((limit - divisor) / step) + 1UL;
|
ulong remainingIterations = ((limit - divisor) / step) + 1UL;
|
||||||
while (true)
|
while (true)
|
||||||
{
|
{
|
||||||
bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0;
|
if (residueStepper.IsAdmissible())
|
||||||
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
|
|
||||||
{
|
{
|
||||||
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
|
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
|
||||||
ulong divisorCycle;
|
ulong divisorCycle;
|
||||||
@@ -792,14 +670,7 @@ private bool CheckDivisors64(
|
|||||||
|
|
||||||
divisor += step;
|
divisor += step;
|
||||||
currentK++;
|
currentK++;
|
||||||
remainder10 = AddMod10(remainder10, step10);
|
residueStepper.Advance();
|
||||||
remainder8 = AddMod8(remainder8, step8);
|
|
||||||
remainder3 = AddMod3(remainder3, step3);
|
|
||||||
remainder7 = AddMod7(remainder7, step7);
|
|
||||||
remainder11 = AddMod11(remainder11, step11);
|
|
||||||
remainder13 = AddMod13(remainder13, step13);
|
|
||||||
remainder17 = AddMod17(remainder17, step17);
|
|
||||||
remainder19 = AddMod19(remainder19, step19);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -870,7 +741,8 @@ private bool CheckDivisors64(
|
|||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
|
|
||||||
BigInteger x = BigInteger.ModPow(baseValue, d, value);
|
var expStepper = new BigIntegerExponentStepper(value, baseValue);
|
||||||
|
BigInteger x = expStepper.Initialize(d);
|
||||||
if (x == 1 || x == value - 1)
|
if (x == 1 || x == value - 1)
|
||||||
{
|
{
|
||||||
continue;
|
continue;
|
||||||
@@ -879,7 +751,8 @@ private bool CheckDivisors64(
|
|||||||
bool witnessFound = true;
|
bool witnessFound = true;
|
||||||
for (int r = 1; r < s; r++)
|
for (int r = 1; r < s; r++)
|
||||||
{
|
{
|
||||||
x = BigInteger.ModPow(x, 2, value);
|
BigInteger targetExponent = d << r;
|
||||||
|
x = expStepper.ComputeNext(targetExponent);
|
||||||
if (x == value - 1)
|
if (x == value - 1)
|
||||||
{
|
{
|
||||||
witnessFound = false;
|
witnessFound = false;
|
||||||
|
|||||||
@@ -2,37 +2,20 @@ using System.Runtime.CompilerServices;
|
|||||||
|
|
||||||
namespace PerfectNumbers.Core;
|
namespace PerfectNumbers.Core;
|
||||||
|
|
||||||
public struct CycleRemainderStepper
|
public struct CycleRemainderStepper(ulong cycleLength)
|
||||||
{
|
{
|
||||||
private readonly ulong _cycleLength;
|
private readonly ulong _cycleLength = cycleLength;
|
||||||
private ulong _previousPrime;
|
private ulong _previousPrime = 0UL;
|
||||||
private ulong _currentRemainder;
|
private ulong _currentRemainder = 0UL;
|
||||||
private bool _hasState;
|
|
||||||
|
|
||||||
public CycleRemainderStepper(ulong cycleLength)
|
|
||||||
{
|
|
||||||
if (cycleLength == 0UL)
|
|
||||||
{
|
|
||||||
throw new ArgumentOutOfRangeException(nameof(cycleLength), "Cycle length must be non-zero.");
|
|
||||||
}
|
|
||||||
|
|
||||||
_cycleLength = cycleLength;
|
|
||||||
_previousPrime = 0UL;
|
|
||||||
_currentRemainder = 0UL;
|
|
||||||
_hasState = false;
|
|
||||||
}
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public void Reset()
|
public void Reset()
|
||||||
{
|
{
|
||||||
// TODO: Inline this reset at the call sites so the hot loops reuse struct reinitialization
|
|
||||||
// measured fastest in MersenneDivisorCycleLengthGpuBenchmarks, avoiding the extra
|
|
||||||
// method call when scanners need to restart stepping.
|
|
||||||
_previousPrime = 0UL;
|
_previousPrime = 0UL;
|
||||||
_currentRemainder = 0UL;
|
_currentRemainder = 0UL;
|
||||||
_hasState = false;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public ulong Initialize(ulong prime)
|
public ulong Initialize(ulong prime)
|
||||||
{
|
{
|
||||||
ulong cycleLength = _cycleLength;
|
ulong cycleLength = _cycleLength;
|
||||||
@@ -44,24 +27,15 @@ public struct CycleRemainderStepper
|
|||||||
// MersenneDivisorCycleLengthGpuBenchmarks winner for both CPU and GPU call sites.
|
// MersenneDivisorCycleLengthGpuBenchmarks winner for both CPU and GPU call sites.
|
||||||
remainder %= cycleLength;
|
remainder %= cycleLength;
|
||||||
}
|
}
|
||||||
|
|
||||||
_previousPrime = prime;
|
_previousPrime = prime;
|
||||||
_currentRemainder = remainder;
|
_currentRemainder = remainder;
|
||||||
_hasState = true;
|
|
||||||
return remainder;
|
return remainder;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public ulong ComputeNext(ulong prime)
|
public ulong ComputeNext(ulong prime)
|
||||||
{
|
{
|
||||||
if (!_hasState)
|
|
||||||
{
|
|
||||||
throw new InvalidOperationException("CycleRemainderStepper must be initialized before computing.");
|
|
||||||
}
|
|
||||||
|
|
||||||
if (prime <= _previousPrime)
|
|
||||||
{
|
|
||||||
throw new ArgumentOutOfRangeException(nameof(prime), "Primes must be processed in strictly increasing order.");
|
|
||||||
}
|
|
||||||
|
|
||||||
ulong cycleLength = _cycleLength;
|
ulong cycleLength = _cycleLength;
|
||||||
ulong delta = prime - _previousPrime;
|
ulong delta = prime - _previousPrime;
|
||||||
_previousPrime = prime;
|
_previousPrime = prime;
|
||||||
@@ -90,13 +64,4 @@ public struct CycleRemainderStepper
|
|||||||
}
|
}
|
||||||
return _currentRemainder;
|
return _currentRemainder;
|
||||||
}
|
}
|
||||||
|
|
||||||
public void CaptureState(out ulong previousPrime, out ulong previousRemainder, out bool hasState)
|
|
||||||
{
|
|
||||||
// TODO: Expose these fields directly once the residue scanners adopt the single-cycle helper
|
|
||||||
// so the hot path can read them without paying for an additional wrapper call per iteration.
|
|
||||||
previousPrime = _previousPrime;
|
|
||||||
previousRemainder = _currentRemainder;
|
|
||||||
hasState = _hasState;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
@@ -1,28 +1,21 @@
|
|||||||
|
using System.Runtime.CompilerServices;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core;
|
namespace PerfectNumbers.Core;
|
||||||
|
|
||||||
/// <summary>
|
/// <summary>
|
||||||
/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path.
|
/// 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.
|
/// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order.
|
||||||
/// </summary>
|
/// </summary>
|
||||||
internal struct ExponentRemainderStepperCpu
|
internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor)
|
||||||
{
|
{
|
||||||
private readonly MontgomeryDivisorData _divisor;
|
private readonly MontgomeryDivisorData _divisor = divisor;
|
||||||
private readonly ulong _modulus;
|
private readonly ulong _modulus = divisor.Modulus;
|
||||||
private readonly ulong _nPrime;
|
private readonly ulong _nPrime = divisor.NPrime;
|
||||||
private readonly ulong _montgomeryOne;
|
private readonly ulong _montgomeryOne = divisor.MontgomeryOne;
|
||||||
public ulong PreviousExponent;
|
private ulong PreviousExponent = 0UL;
|
||||||
private ulong _currentMontgomery;
|
private ulong _currentMontgomery = divisor.MontgomeryOne;
|
||||||
|
|
||||||
public ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor)
|
|
||||||
{
|
|
||||||
_divisor = divisor;
|
|
||||||
_modulus = divisor.Modulus;
|
|
||||||
_nPrime = divisor.NPrime;
|
|
||||||
_montgomeryOne = divisor.MontgomeryOne;
|
|
||||||
PreviousExponent = 0UL;
|
|
||||||
_currentMontgomery = divisor.MontgomeryOne;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor)
|
public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor)
|
||||||
=> _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne;
|
=> _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne;
|
||||||
|
|
||||||
@@ -30,6 +23,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// Clears the cached Montgomery residue so the stepper can be reused for another divisor.
|
/// 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.
|
/// Callers must reinitialize the stepper before consuming any residues after a reset.
|
||||||
/// </summary>
|
/// </summary>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public void Reset()
|
public void Reset()
|
||||||
{
|
{
|
||||||
PreviousExponent = 0UL;
|
PreviousExponent = 0UL;
|
||||||
@@ -42,6 +36,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// </summary>
|
/// </summary>
|
||||||
/// <param name="exponent">First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata.</param>
|
/// <param name="exponent">First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata.</param>
|
||||||
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
|
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public ulong InitializeCpu(ulong exponent)
|
public ulong InitializeCpu(ulong exponent)
|
||||||
{
|
{
|
||||||
InitializeCpuState(exponent);
|
InitializeCpuState(exponent);
|
||||||
@@ -53,6 +48,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// </summary>
|
/// </summary>
|
||||||
/// <param name="exponent">First exponent evaluated for the current divisor. Must match the divisor passed to the constructor.</param>
|
/// <param name="exponent">First exponent evaluated for the current divisor. Must match the divisor passed to the constructor.</param>
|
||||||
/// <returns><see langword="true"/> when the seed exponent produces a Montgomery residue equal to one.</returns>
|
/// <returns><see langword="true"/> when the seed exponent produces a Montgomery residue equal to one.</returns>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public bool InitializeCpuIsUnity(ulong exponent)
|
public bool InitializeCpuIsUnity(ulong exponent)
|
||||||
{
|
{
|
||||||
InitializeCpuState(exponent);
|
InitializeCpuState(exponent);
|
||||||
@@ -65,6 +61,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// </summary>
|
/// </summary>
|
||||||
/// <param name="exponent">Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent.</param>
|
/// <param name="exponent">Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent.</param>
|
||||||
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
|
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public ulong ComputeNextCpu(ulong exponent)
|
public ulong ComputeNextCpu(ulong exponent)
|
||||||
{
|
{
|
||||||
ulong delta = exponent - PreviousExponent;
|
ulong delta = exponent - PreviousExponent;
|
||||||
@@ -84,6 +81,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// </summary>
|
/// </summary>
|
||||||
/// <param name="exponent">Next exponent in the ascending sequence.</param>
|
/// <param name="exponent">Next exponent in the ascending sequence.</param>
|
||||||
/// <returns><see langword="true"/> when the canonical residue equals one.</returns>
|
/// <returns><see langword="true"/> when the canonical residue equals one.</returns>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public bool ComputeNextIsUnity(ulong exponent)
|
public bool ComputeNextIsUnity(ulong exponent)
|
||||||
{
|
{
|
||||||
ulong delta = exponent - PreviousExponent;
|
ulong delta = exponent - PreviousExponent;
|
||||||
@@ -103,6 +101,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// <param name="montgomeryResult">Residue in Montgomery form produced externally for the configured divisor.</param>
|
/// <param name="montgomeryResult">Residue in Montgomery form produced externally for the configured divisor.</param>
|
||||||
/// <param name="isUnity">Outputs whether the supplied residue equals Montgomery one.</param>
|
/// <param name="isUnity">Outputs whether the supplied residue equals Montgomery one.</param>
|
||||||
/// <returns><see langword="true"/> when the residue is accepted and cached.</returns>
|
/// <returns><see langword="true"/> when the residue is accepted and cached.</returns>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity)
|
public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity)
|
||||||
{
|
{
|
||||||
_currentMontgomery = montgomeryResult;
|
_currentMontgomery = montgomeryResult;
|
||||||
@@ -119,6 +118,7 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
/// <param name="montgomeryDelta">Montgomery residue delta corresponding to the exponent gap.</param>
|
/// <param name="montgomeryDelta">Montgomery residue delta corresponding to the exponent gap.</param>
|
||||||
/// <param name="isUnity">Outputs whether the resulting residue equals Montgomery one.</param>
|
/// <param name="isUnity">Outputs whether the resulting residue equals Montgomery one.</param>
|
||||||
/// <returns><see langword="true"/> when the advance succeeds.</returns>
|
/// <returns><see langword="true"/> when the advance succeeds.</returns>
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity)
|
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
|
// TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead
|
||||||
@@ -130,11 +130,13 @@ internal struct ExponentRemainderStepperCpu
|
|||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
private void InitializeCpuState(ulong exponent)
|
private void InitializeCpuState(ulong exponent)
|
||||||
{
|
{
|
||||||
_currentMontgomery = exponent.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
|
_currentMontgomery = exponent.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor);
|
||||||
PreviousExponent = exponent;
|
PreviousExponent = exponent;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
private readonly ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime);
|
private readonly ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime);
|
||||||
}
|
}
|
||||||
@@ -11,7 +11,7 @@ internal struct ExponentRemainderStepperGpu
|
|||||||
{
|
{
|
||||||
private readonly ulong _modulus;
|
private readonly ulong _modulus;
|
||||||
private readonly ReadOnlyGpuUInt128 _modulusWide;
|
private readonly ReadOnlyGpuUInt128 _modulusWide;
|
||||||
public ulong PreviousExponent;
|
private ulong PreviousExponent;
|
||||||
private GpuUInt128 _currentResidue;
|
private GpuUInt128 _currentResidue;
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
|||||||
@@ -5,6 +5,7 @@ using System.Runtime.CompilerServices;
|
|||||||
using System.Runtime.InteropServices;
|
using System.Runtime.InteropServices;
|
||||||
using ILGPU;
|
using ILGPU;
|
||||||
using ILGPU.Runtime;
|
using ILGPU.Runtime;
|
||||||
|
using PerfectNumbers.Core;
|
||||||
using PerfectNumbers.Core.Gpu.Accelerators;
|
using PerfectNumbers.Core.Gpu.Accelerators;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core.Gpu;
|
namespace PerfectNumbers.Core.Gpu;
|
||||||
@@ -253,16 +254,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
|
ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
|
||||||
ulong startK = minK < 1UL ? 1UL : minK;
|
ulong startK = minK < 1UL ? 1UL : minK;
|
||||||
|
|
||||||
ulong stepValue = (ulong)twoP128;
|
|
||||||
byte step10 = (byte)(stepValue % 10UL);
|
|
||||||
byte step8 = (byte)(stepValue & 7UL);
|
|
||||||
byte step3 = (byte)(stepValue % 3UL);
|
|
||||||
byte step7 = (byte)(stepValue % 7UL);
|
|
||||||
byte step11 = (byte)(stepValue % 11UL);
|
|
||||||
|
|
||||||
LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
|
|
||||||
ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);
|
|
||||||
|
|
||||||
Span<ulong> divisorStorage = divisors.AsSpan();
|
Span<ulong> divisorStorage = divisors.AsSpan();
|
||||||
Span<ulong> exponentStorage = exponents.AsSpan();
|
Span<ulong> exponentStorage = exponents.AsSpan();
|
||||||
Span<GpuDivisorPartialData> divisorDataStorage = divisorData.AsSpan();
|
Span<GpuDivisorPartialData> divisorDataStorage = divisorData.AsSpan();
|
||||||
@@ -289,12 +280,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
ulong prime,
|
ulong prime,
|
||||||
UInt128 twoP128,
|
UInt128 twoP128,
|
||||||
UInt128 allowedMax128,
|
UInt128 allowedMax128,
|
||||||
byte step10,
|
|
||||||
byte step8,
|
|
||||||
byte step3,
|
|
||||||
byte step7,
|
|
||||||
byte step11,
|
|
||||||
ushort decimalMask,
|
|
||||||
ArrayView1D<GpuDivisorPartialData, Stride1D.Dense> divisorDataView,
|
ArrayView1D<GpuDivisorPartialData, Stride1D.Dense> divisorDataView,
|
||||||
ArrayView1D<int, Stride1D.Dense> offsetView,
|
ArrayView1D<int, Stride1D.Dense> offsetView,
|
||||||
ArrayView1D<int, Stride1D.Dense> countView,
|
ArrayView1D<int, Stride1D.Dense> countView,
|
||||||
@@ -329,6 +314,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)twoP128, (GpuUInt128)startDivisor128);
|
||||||
UInt128 rangeLimit128 = (twoP128 * rangeEndK) + UInt128.One;
|
UInt128 rangeLimit128 = (twoP128 * rangeEndK) + UInt128.One;
|
||||||
if (rangeLimit128 > allowedMax128)
|
if (rangeLimit128 > allowedMax128)
|
||||||
{
|
{
|
||||||
@@ -338,12 +324,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
UInt128 rangeCount128 = ((rangeLimit128 - startDivisor128) / twoP128) + UInt128.One;
|
UInt128 rangeCount128 = ((rangeLimit128 - startDivisor128) / twoP128) + UInt128.One;
|
||||||
ulong rangeRemaining = rangeCount128 > ulong.MaxValue ? ulong.MaxValue : (ulong)rangeCount128;
|
ulong rangeRemaining = rangeCount128 > ulong.MaxValue ? ulong.MaxValue : (ulong)rangeCount128;
|
||||||
UInt128 currentRangeDivisor128 = startDivisor128;
|
UInt128 currentRangeDivisor128 = startDivisor128;
|
||||||
ulong currentRangeDivisorValue = (ulong)currentRangeDivisor128;
|
|
||||||
byte rangeRemainder10 = (byte)(currentRangeDivisorValue % 10UL);
|
|
||||||
byte rangeRemainder8 = (byte)(currentRangeDivisorValue & 7UL);
|
|
||||||
byte rangeRemainder3 = (byte)(currentRangeDivisorValue % 3UL);
|
|
||||||
byte rangeRemainder7 = (byte)(currentRangeDivisorValue % 7UL);
|
|
||||||
byte rangeRemainder11 = (byte)(currentRangeDivisorValue % 11UL);
|
|
||||||
|
|
||||||
while (rangeRemaining > 0UL && !composite)
|
while (rangeRemaining > 0UL && !composite)
|
||||||
{
|
{
|
||||||
@@ -356,16 +336,11 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
UInt128 nextDivisor128 = currentRangeDivisor128;
|
UInt128 nextDivisor128 = currentRangeDivisor128;
|
||||||
int admissibleCount = 0;
|
int admissibleCount = 0;
|
||||||
|
|
||||||
byte localRemainder10 = rangeRemainder10;
|
var localStepper = residueStepper;
|
||||||
byte localRemainder8 = rangeRemainder8;
|
|
||||||
byte localRemainder3 = rangeRemainder3;
|
|
||||||
byte localRemainder7 = rangeRemainder7;
|
|
||||||
byte localRemainder11 = rangeRemainder11;
|
|
||||||
|
|
||||||
for (int i = 0; i < chunkCount; i++)
|
for (int i = 0; i < chunkCount; i++)
|
||||||
{
|
{
|
||||||
bool passesSmallModuli = localRemainder3 != 0 && localRemainder7 != 0 && localRemainder11 != 0;
|
if (localStepper.IsAdmissible())
|
||||||
if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0)
|
|
||||||
{
|
{
|
||||||
ulong divisorValue = (ulong)nextDivisor128;
|
ulong divisorValue = (ulong)nextDivisor128;
|
||||||
MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue);
|
MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue);
|
||||||
@@ -389,18 +364,10 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
}
|
}
|
||||||
|
|
||||||
nextDivisor128 += twoP128;
|
nextDivisor128 += twoP128;
|
||||||
localRemainder10 = AddMod10(localRemainder10, step10);
|
localStepper.Advance();
|
||||||
localRemainder8 = AddMod8(localRemainder8, step8);
|
|
||||||
localRemainder3 = AddMod3(localRemainder3, step3);
|
|
||||||
localRemainder7 = AddMod7(localRemainder7, step7);
|
|
||||||
localRemainder11 = AddMod11(localRemainder11, step11);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
rangeRemainder10 = localRemainder10;
|
residueStepper = localStepper;
|
||||||
rangeRemainder8 = localRemainder8;
|
|
||||||
rangeRemainder3 = localRemainder3;
|
|
||||||
rangeRemainder7 = localRemainder7;
|
|
||||||
rangeRemainder11 = localRemainder11;
|
|
||||||
|
|
||||||
processedCountLocal += (ulong)chunkCount;
|
processedCountLocal += (ulong)chunkCount;
|
||||||
|
|
||||||
@@ -456,12 +423,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
prime,
|
prime,
|
||||||
twoP128,
|
twoP128,
|
||||||
allowedMax128,
|
allowedMax128,
|
||||||
step10,
|
|
||||||
step8,
|
|
||||||
step3,
|
|
||||||
step7,
|
|
||||||
step11,
|
|
||||||
decimalMask,
|
|
||||||
divisorDataView,
|
divisorDataView,
|
||||||
offsetView,
|
offsetView,
|
||||||
countView,
|
countView,
|
||||||
@@ -502,12 +463,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
prime,
|
prime,
|
||||||
twoP128,
|
twoP128,
|
||||||
allowedMax128,
|
allowedMax128,
|
||||||
step10,
|
|
||||||
step8,
|
|
||||||
step3,
|
|
||||||
step7,
|
|
||||||
step11,
|
|
||||||
decimalMask,
|
|
||||||
divisorDataView,
|
divisorDataView,
|
||||||
offsetView,
|
offsetView,
|
||||||
countView,
|
countView,
|
||||||
@@ -550,69 +505,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
|
|||||||
return computedCycle;
|
return computedCycle;
|
||||||
}
|
}
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
||||||
private static byte AddMod3(byte value, byte delta)
|
|
||||||
{
|
|
||||||
const int Modulus = 3;
|
|
||||||
int sum = value + delta;
|
|
||||||
|
|
||||||
if (sum >= Modulus)
|
|
||||||
{
|
|
||||||
sum -= Modulus;
|
|
||||||
}
|
|
||||||
|
|
||||||
return (byte)sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
||||||
private static byte AddMod7(byte value, byte delta)
|
|
||||||
{
|
|
||||||
const int Modulus = 7;
|
|
||||||
int sum = value + delta;
|
|
||||||
|
|
||||||
if (sum >= Modulus)
|
|
||||||
{
|
|
||||||
sum -= Modulus;
|
|
||||||
}
|
|
||||||
|
|
||||||
return (byte)sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
||||||
private static byte AddMod10(byte value, byte delta)
|
|
||||||
{
|
|
||||||
const int Modulus = 10;
|
|
||||||
int sum = value + delta;
|
|
||||||
|
|
||||||
if (sum >= Modulus)
|
|
||||||
{
|
|
||||||
sum -= Modulus;
|
|
||||||
}
|
|
||||||
|
|
||||||
return (byte)sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
||||||
private static byte AddMod11(byte value, byte delta)
|
|
||||||
{
|
|
||||||
const int Modulus = 11;
|
|
||||||
int sum = value + delta;
|
|
||||||
|
|
||||||
if (sum >= Modulus)
|
|
||||||
{
|
|
||||||
sum -= Modulus;
|
|
||||||
}
|
|
||||||
|
|
||||||
return (byte)sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
||||||
private static byte AddMod8(byte value, byte delta)
|
|
||||||
{
|
|
||||||
return (byte)((value + delta) & 7);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime)
|
private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime)
|
||||||
{
|
{
|
||||||
|
|||||||
@@ -3,50 +3,32 @@ using PerfectNumbers.Core.Gpu;
|
|||||||
|
|
||||||
namespace PerfectNumbers.Core;
|
namespace PerfectNumbers.Core;
|
||||||
|
|
||||||
internal struct MersenneDivisorResidueStepper
|
[method: MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
internal struct MersenneDivisorResidueStepper(ulong prime, in GpuUInt128 step, in GpuUInt128 firstDivisor)
|
||||||
{
|
{
|
||||||
private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
|
private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
|
||||||
private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
|
private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
|
||||||
|
|
||||||
private readonly ushort decimalMask;
|
private readonly ushort decimalMask = (prime & 3UL) == 3UL ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
|
||||||
private readonly byte step10;
|
private readonly byte step10 = ComputeModulo(step, 10, multiplier: 6);
|
||||||
private readonly byte step8;
|
private readonly byte step8 = (byte)(step.Low % 8UL);
|
||||||
private readonly byte step5;
|
private readonly byte step5 = ComputeModulo(step, 5, multiplier: 1);
|
||||||
private readonly byte step3;
|
private readonly byte step3 = ComputeModulo(step, 3, multiplier: 1);
|
||||||
private readonly byte step7;
|
private readonly byte step7 = ComputeModulo(step, 7, multiplier: 2);
|
||||||
private readonly byte step11;
|
private readonly byte step11 = ComputeModulo(step, 11, multiplier: 5);
|
||||||
|
private readonly byte step13 = ComputeModulo(step, 13, multiplier: 3);
|
||||||
|
private readonly byte step17 = ComputeModulo(step, 17, multiplier: 1);
|
||||||
|
private readonly byte step19 = ComputeModulo(step, 19, multiplier: 17);
|
||||||
|
|
||||||
public byte Remainder10 { get; private set; }
|
private byte Remainder10 = ComputeModulo(firstDivisor, 10, multiplier: 6);
|
||||||
|
private byte Remainder8 = (byte)(firstDivisor.Low % 8UL);
|
||||||
public byte Remainder8 { get; private set; }
|
private byte Remainder5 = ComputeModulo(firstDivisor, 5, multiplier: 1);
|
||||||
|
private byte Remainder3 = ComputeModulo(firstDivisor, 3, multiplier: 1);
|
||||||
public byte Remainder5 { get; private set; }
|
private byte Remainder7 = ComputeModulo(firstDivisor, 7, multiplier: 2);
|
||||||
|
private byte Remainder11 = ComputeModulo(firstDivisor, 11, multiplier: 5);
|
||||||
public byte Remainder3 { get; private set; }
|
private byte Remainder13 = ComputeModulo(firstDivisor, 13, multiplier: 3);
|
||||||
|
private byte Remainder17 = ComputeModulo(firstDivisor, 17, multiplier: 1);
|
||||||
public byte Remainder7 { get; private set; }
|
private byte Remainder19 = ComputeModulo(firstDivisor, 19, multiplier: 17);
|
||||||
|
|
||||||
public byte Remainder11 { get; private set; }
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
||||||
public MersenneDivisorResidueStepper(ulong prime, in GpuUInt128 step, in GpuUInt128 firstDivisor)
|
|
||||||
{
|
|
||||||
decimalMask = (prime & 3UL) == 3UL ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
|
|
||||||
|
|
||||||
step10 = ComputeModulo(step, 10, multiplier: 6);
|
|
||||||
step8 = (byte)(step.Low % 8UL);
|
|
||||||
step5 = ComputeModulo(step, 5, multiplier: 1);
|
|
||||||
step3 = ComputeModulo(step, 3, multiplier: 1);
|
|
||||||
step7 = ComputeModulo(step, 7, multiplier: 2);
|
|
||||||
step11 = ComputeModulo(step, 11, multiplier: 5);
|
|
||||||
|
|
||||||
Remainder10 = ComputeModulo(firstDivisor, 10, multiplier: 6);
|
|
||||||
Remainder8 = (byte)(firstDivisor.Low % 8UL);
|
|
||||||
Remainder5 = ComputeModulo(firstDivisor, 5, multiplier: 1);
|
|
||||||
Remainder3 = ComputeModulo(firstDivisor, 3, multiplier: 1);
|
|
||||||
Remainder7 = ComputeModulo(firstDivisor, 7, multiplier: 2);
|
|
||||||
Remainder11 = ComputeModulo(firstDivisor, 11, multiplier: 5);
|
|
||||||
}
|
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public readonly bool IsAdmissible()
|
public readonly bool IsAdmissible()
|
||||||
@@ -61,7 +43,7 @@ internal struct MersenneDivisorResidueStepper
|
|||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
if (Remainder3 == 0 || Remainder5 == 0 || Remainder7 == 0 || Remainder11 == 0)
|
if (Remainder3 == 0 || Remainder5 == 0 || Remainder7 == 0 || Remainder11 == 0 || Remainder13 == 0 || Remainder17 == 0 || Remainder19 == 0)
|
||||||
{
|
{
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
@@ -78,6 +60,9 @@ internal struct MersenneDivisorResidueStepper
|
|||||||
Remainder3 = AddMod(Remainder3, step3, 3);
|
Remainder3 = AddMod(Remainder3, step3, 3);
|
||||||
Remainder7 = AddMod(Remainder7, step7, 7);
|
Remainder7 = AddMod(Remainder7, step7, 7);
|
||||||
Remainder11 = AddMod(Remainder11, step11, 11);
|
Remainder11 = AddMod(Remainder11, step11, 11);
|
||||||
|
Remainder13 = AddMod(Remainder13, step13, 13);
|
||||||
|
Remainder17 = AddMod(Remainder17, step17, 17);
|
||||||
|
Remainder19 = AddMod(Remainder19, step19, 19);
|
||||||
}
|
}
|
||||||
|
|
||||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
|||||||
@@ -1,19 +1,17 @@
|
|||||||
|
using System.Runtime.CompilerServices;
|
||||||
|
|
||||||
namespace PerfectNumbers.Core;
|
namespace PerfectNumbers.Core;
|
||||||
|
|
||||||
public class MersenneResidueStepper
|
public struct MersenneResidueStepper
|
||||||
{
|
{
|
||||||
private UInt128 m;
|
private readonly UInt128 m;
|
||||||
private ulong p;
|
private ulong p;
|
||||||
private UInt128 pow2_mod;
|
private UInt128 pow2_mod;
|
||||||
private UInt128 M_mod;
|
private UInt128 M_mod;
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public MersenneResidueStepper(UInt128 modulus, ulong p0, UInt128 M_mod_m_at_p0)
|
public MersenneResidueStepper(UInt128 modulus, ulong p0, UInt128 M_mod_m_at_p0)
|
||||||
{
|
{
|
||||||
if (modulus == 0)
|
|
||||||
{
|
|
||||||
throw new ArgumentException("Modulus must be positive.");
|
|
||||||
}
|
|
||||||
|
|
||||||
m = modulus;
|
m = modulus;
|
||||||
p = p0;
|
p = p0;
|
||||||
|
|
||||||
@@ -32,8 +30,10 @@ public class MersenneResidueStepper
|
|||||||
|
|
||||||
// TODO: Inline this accessor into callers so the hot residue path reads the backing fields directly
|
// TODO: Inline this accessor into callers so the hot residue path reads the backing fields directly
|
||||||
// instead of jumping through a wrapper method every iteration.
|
// instead of jumping through a wrapper method every iteration.
|
||||||
public (ulong, UInt128) State() => (p, M_mod);
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
|
public readonly (ulong, UInt128) State() => (p, M_mod);
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public (ulong, UInt128) Step(long delta)
|
public (ulong, UInt128) Step(long delta)
|
||||||
{
|
{
|
||||||
if (delta == 0)
|
if (delta == 0)
|
||||||
@@ -67,6 +67,7 @@ public class MersenneResidueStepper
|
|||||||
return State();
|
return State();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
public (ulong, UInt128) To(ulong p_target)
|
public (ulong, UInt128) To(ulong p_target)
|
||||||
{
|
{
|
||||||
// TODO: Collapse this wrapper once callers can invoke Step directly; the extra indirection shows up
|
// TODO: Collapse this wrapper once callers can invoke Step directly; the extra indirection shows up
|
||||||
@@ -76,6 +77,7 @@ public class MersenneResidueStepper
|
|||||||
return Step(delta);
|
return Step(delta);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||||
private static UInt128 ModInverse(UInt128 a, UInt128 m)
|
private static UInt128 ModInverse(UInt128 a, UInt128 m)
|
||||||
{
|
{
|
||||||
UInt128 m0 = m;
|
UInt128 m0 = m;
|
||||||
|
|||||||
Reference in new issue
Block a user