mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

PrimeOrderCalculatorAccelerator.Rent is now created outside of the loops and re-used for all calculations.

Browse Source
This commit is contained in:
mlemanczyk committed 2025-11-24 22:51:53 +01:00
1 parent 57a50e2a37
commit b9ca95bf8c
53 files changed
+4639 -4229
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs
+6 -3
View File
@@ -1,11 +1,14 @@
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Cpu;
internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelerator gpu) : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{
public void Reset()
public void Reset()
{
}
// TODO: Convert this to static method
public void CheckDivisor(
ulong divisor,
in MontgomeryDivisorData divisorData,
@@ -29,7 +32,7 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi
if (divisorCycle == 0UL)
{
divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
divisorCycle = DivisorCycleCache.Shared.GetCycleLength(gpu, divisor);
if (divisorCycle == 0UL)
{
// DivisorCycleCache guarantees a positive cycle for divisors greater than one.
PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+11 -14
View File
@@ -4,7 +4,7 @@ using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Cpu;
public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester
public sealed class MersenneNumberDivisorByDivisorCpuTester() : IMersenneNumberDivisorByDivisorTester
{
private ulong _divisorLimit;
private int _batchSize = 1_024;
@@ -42,7 +42,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public ulong GetAllowedMaxDivisor(ulong prime) => ComputeAllowedMaxDivisor(prime, _divisorLimit);
public bool IsPrime(ulong prime, out bool divisorsExhausted, out ulong divisor)
public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor)
{
ulong allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
@@ -58,6 +58,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
bool processedAll;
bool composite = CheckDivisors(
gpu,
prime,
allowedMax,
out processedAll,
@@ -84,7 +85,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
}
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession()
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu)
{
MersenneCpuDivisorScanSession? session = ThreadStaticPools.RentMersenneCpuDivisorSession();
if (session is not null)
@@ -92,16 +93,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return session;
}
return new MersenneCpuDivisorScanSession();
return new MersenneCpuDivisorScanSession(gpu);
}
private bool CheckDivisors(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong allowedMax,
out bool processedAll,
out ulong foundDivisor)
{
foundDivisor = 0UL;
// The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs.
// Keeping the guard commented out documents the reasoning for benchmarks and tests.
// if (allowedMax < 3UL)
@@ -171,6 +172,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
remainder11 = (byte)(divisorLow % 11UL);
return CheckDivisors64(
gpu,
prime,
stepLow,
limit.Low,
@@ -203,7 +205,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
var divisorPool = MontgomeryDivisorDataPool.Shared;
while (divisor.CompareTo(limit) <= 0)
{
@@ -225,6 +226,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
// Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles
// for the exponents scanned here, so skip the unused cache fallback and compute directly.
divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
gpu,
candidate,
divisorData,
skipPrimeOrderHeuristic: primeOrderFailed);
@@ -241,7 +243,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
// the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p.
foundDivisor = candidate;
processedAll = true;
PrimeOrderCalculatorAccelerator.Return(gpu);
return true;
}
@@ -259,13 +260,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
remainder11 = AddMod11(remainder11, step11);
}
PrimeOrderCalculatorAccelerator.Return(gpu);
processedAll = true;
foundDivisor = 0UL;
return false;
}
private static bool CheckDivisors64(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong step,
ulong limit,
@@ -285,7 +286,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
out ulong foundDivisor)
{
bool canAdvance = step <= limit;
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
if (!canAdvance)
{
@@ -304,6 +304,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
out bool primeOrderFailed) || computedCycle == 0UL)
{
divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
gpu,
divisor,
divisorData,
skipPrimeOrderHeuristic: primeOrderFailed);
@@ -319,7 +320,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
{
foundDivisor = divisor;
processedAll = true;
PrimeOrderCalculatorAccelerator.Return(gpu);
return true;
}
@@ -332,7 +332,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
processedAll = true;
foundDivisor = 0UL;
PrimeOrderCalculatorAccelerator.Return(gpu);
return false;
}
@@ -340,7 +339,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
{
processedAll = true;
foundDivisor = 0UL;
PrimeOrderCalculatorAccelerator.Return(gpu);
return false;
}
@@ -361,6 +359,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
out bool primeOrderFailed) || computedCycle == 0UL)
{
divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
gpu,
divisor,
divisorData,
skipPrimeOrderHeuristic: primeOrderFailed);
@@ -375,7 +374,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
{
foundDivisor = divisor;
processedAll = true;
PrimeOrderCalculatorAccelerator.Return(gpu);
return true;
}
@@ -390,7 +388,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
{
processedAll = true;
foundDivisor = 0UL;
PrimeOrderCalculatorAccelerator.Return(gpu);
return false;
}
PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs
+80 -78
View File
@@ -1,90 +1,92 @@
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Cpu;
public class MersenneNumberIncrementalCpuTester(GpuKernelType kernelType)
{
private static readonly MersenneDivisorCycles _cycles = MersenneDivisorCycles.Shared;
private readonly GpuKernelType _kernelType = kernelType;
private static readonly MersenneDivisorCycles _cycles = MersenneDivisorCycles.Shared;
private readonly GpuKernelType _kernelType = kernelType;
public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
{
UInt128 k = UInt128.One;
UInt128 one = k;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
// Use residue automaton always on CPU
var auto = new MersenneResidueAutomaton(exponent);
public void Scan(PrimeOrderCalculatorAccelerator gpu, ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
{
UInt128 k = UInt128.One;
UInt128 one = k;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
// Use residue automaton always on CPU
var auto = new MersenneResidueAutomaton(exponent);
bool reject;
bool shouldCheck;
UInt128 q = auto.CurrentQ();
UInt128 divPow;
UInt128 halfPow;
UInt128 phi;
UInt128 qCycle;
bool reject;
bool shouldCheck;
UInt128 q = auto.CurrentQ();
UInt128 divPow;
UInt128 halfPow;
UInt128 phi;
UInt128 qCycle;
ulong div;
ulong mod8;
ulong phi64;
ulong div;
ulong mod8;
ulong phi64;
while (k <= maxK && Volatile.Read(ref isPrime))
{
// TODO: Replace this direct GetCycle call with DivisorCycleCache.Lookup so we reuse the single snapshot block
// and compute missing cycles on the configured device without queuing additional block generation.
qCycle = MersenneDivisorCycles.GetCycle(q);
shouldCheck = auto.Mod10R != 5UL;
if (shouldCheck)
{
mod8 = auto.Mod8R;
if (mod8 != 1UL && mod8 != 7UL)
{
shouldCheck = false;
}
else if (auto.Mod3R == 0UL || auto.Mod5R == 0UL)
{
// TODO: Swap these Mod3/Mod5 zero checks to the cached residue tables once the automaton exposes
// the benchmarked lookup-based helpers so CPU scans avoid runtime modulo instructions.
shouldCheck = false;
}
}
while (k <= maxK && Volatile.Read(ref isPrime))
{
// TODO: Replace this direct GetCycle call with DivisorCycleCache.Lookup so we reuse the single snapshot block
// and compute missing cycles on the configured device without queuing additional block generation.
qCycle = MersenneDivisorCycles.GetCycle(gpu, q);
shouldCheck = auto.Mod10R != 5UL;
if (shouldCheck)
{
mod8 = auto.Mod8R;
if (mod8 != 1UL && mod8 != 7UL)
{
shouldCheck = false;
}
else if (auto.Mod3R == 0UL || auto.Mod5R == 0UL)
{
// TODO: Swap these Mod3/Mod5 zero checks to the cached residue tables once the automaton exposes
// the benchmarked lookup-based helpers so CPU scans avoid runtime modulo instructions.
shouldCheck = false;
}
}
if (shouldCheck)
{
if (_kernelType == GpuKernelType.Pow2Mod)
{
reject = exponent.PowModWithCycle(q, qCycle) == UInt128.One;
}
else
{
phi = q - one;
reject = false;
if (phi <= ulong.MaxValue)
{
phi64 = (ulong)phi;
if (phi64.PowModWithCycle(q, qCycle) == one)
{
halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - one;
if (halfPow.BinaryGcd(q) == one)
{
div = phi64.FastDiv64(exponent, divMul);
divPow = div.PowModWithCycle(q, qCycle) - one;
if (divPow.BinaryGcd(q) == one)
{
reject = true;
}
}
}
}
}
if (shouldCheck)
{
if (_kernelType == GpuKernelType.Pow2Mod)
{
reject = exponent.PowModWithCycle(q, qCycle) == UInt128.One;
}
else
{
phi = q - one;
reject = false;
if (phi <= ulong.MaxValue)
{
phi64 = (ulong)phi;
if (phi64.PowModWithCycle(q, qCycle) == one)
{
halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - one;
if (halfPow.BinaryGcd(q) == one)
{
div = phi64.FastDiv64(exponent, divMul);
divPow = div.PowModWithCycle(q, qCycle) - one;
if (divPow.BinaryGcd(q) == one)
{
reject = true;
}
}
}
}
}
if (reject && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
break;
}
}
if (reject && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
break;
}
}
k += 1UL;
auto.Next();
q = auto.CurrentQ();
}
}
k += 1UL;
auto.Next();
q = auto.CurrentQ();
}
}
}
PerfectNumbers.Core/Cpu/MersenneNumberLucasLehmerCpuTester.cs
+37 -35
View File
@@ -1,44 +1,46 @@
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core;
public class MersenneNumberLucasLehmerCpuTester
{
public bool IsPrime(ulong exponent)
{
// Early rejections aligned with incremental/order sieves, but safe for small p:
// - If 3 | p and p != 3, then 7 | M_p -> composite.
// - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast.
// - If p is divisible by divisors specific to numbers ending with 1 or 7, reject fast.
// TODO: Replace these `%` checks with Mod3/Mod5/Mod7/Mod11 helpers once Lucas–Lehmer CPU filtering
// shares the benchmarked bitmask implementations and avoids slow modulo instructions in the hot path.
if (
((exponent % 3UL) == 0UL && exponent != 3UL) ||
((exponent % 5UL) == 0UL && exponent != 5UL) ||
((exponent % 7UL) == 0UL && exponent != 7UL) ||
((exponent % 11UL) == 0UL && exponent != 11UL)
)
{
return false;
}
public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong exponent)
{
// Early rejections aligned with incremental/order sieves, but safe for small p:
// - If 3 | p and p != 3, then 7 | M_p -> composite.
// - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast.
// - If p is divisible by divisors specific to numbers ending with 1 or 7, reject fast.
// TODO: Replace these `%` checks with Mod3/Mod5/Mod7/Mod11 helpers once Lucas–Lehmer CPU filtering
// shares the benchmarked bitmask implementations and avoids slow modulo instructions in the hot path.
if (
((exponent % 3UL) == 0UL && exponent != 3UL) ||
((exponent % 5UL) == 0UL && exponent != 5UL) ||
((exponent % 7UL) == 0UL && exponent != 7UL) ||
((exponent % 11UL) == 0UL && exponent != 11UL)
)
{
return false;
}
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
{
return false;
}
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne(gpu))
{
return false;
}
UInt128 s = UInt128Numbers.Four, two = UInt128Numbers.Two,
m = (UInt128.One << (int)exponent) - UInt128.One;
UInt128 s = UInt128Numbers.Four, two = UInt128Numbers.Two,
m = (UInt128.One << (int)exponent) - UInt128.One;
UInt128 limit = exponent - UInt128Numbers.Two;
for (UInt128 i = UInt128.Zero; i < limit; i++)
{
// Lucas–Lehmer squares the evolving residue, so we stay on the Montgomery powmod path rather than the pow2-specific helper.
s = s.PowModWithCycle(two, MersenneDivisorCycles.GetCycle(m)) - two;
if (s < UInt128.Zero)
{
s += m;
}
}
UInt128 limit = exponent - UInt128Numbers.Two;
for (UInt128 i = UInt128.Zero; i < limit; i++)
{
// Lucas–Lehmer squares the evolving residue, so we stay on the Montgomery powmod path rather than the pow2-specific helper.
s = s.PowModWithCycle(two, MersenneDivisorCycles.GetCycle(gpu, m)) - two;
if (s < UInt128.Zero)
{
s += m;
}
}
return s == UInt128.Zero;
}
return s == UInt128.Zero;
}
}
PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs
+35 -33
View File
@@ -1,10 +1,12 @@
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Cpu;
public class MersenneNumberOrderCpuTester(GpuKernelType kernelType)
{
private readonly GpuKernelType _kernelType = kernelType;
public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
public void Scan(PrimeOrderCalculatorAccelerator gpu, ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
{
UInt128 k = 1UL;
var auto = new MersenneResidueAutomaton(exponent);
@@ -13,10 +15,10 @@ public class MersenneNumberOrderCpuTester(GpuKernelType kernelType)
while (k <= maxK && Volatile.Read(ref isPrime))
{
qCycle = MersenneDivisorCycles.GetCycle(q);
// TODO: When this lookup misses the snapshot, invoke the configured device to compute the
// single required cycle on demand without persisting it so the CPU order path retains the
// cycle-stepping speedups while honoring the no-extra-cache constraint for large divisors.
qCycle = MersenneDivisorCycles.GetCycle(gpu, q);
// TODO: When this lookup misses the snapshot, invoke the configured device to compute the
// single required cycle on demand without persisting it so the CPU order path retains the
// cycle-stepping speedups while honoring the no-extra-cache constraint for large divisors.
int allowMask2 = lastDigit == LastDigit.Seven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10;
bool shouldCheck = ((allowMask2 >> (int)auto.Mod10R) & 1) != 0;
if (shouldCheck)
@@ -33,33 +35,33 @@ public class MersenneNumberOrderCpuTester(GpuKernelType kernelType)
}
if (shouldCheck)
{
if (_kernelType == GpuKernelType.Pow2Mod)
{
if (exponent.PowModWithCycle(q, qCycle) == 1UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
else
{
UInt128 phi = q - 1UL;
if (phi <= ulong.MaxValue)
{
ulong phi64 = (ulong)phi;
if (phi64.PowModWithCycle(q, qCycle) == 1UL)
{
UInt128 halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - 1UL;
if (halfPow.BinaryGcd(q) == 1UL)
{
ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL;
ulong div = phi64.FastDiv64(exponent, divMul);
UInt128 divPow = div.PowModWithCycle(q, qCycle) - 1UL;
if (divPow.BinaryGcd(q) == 1UL)
{
Volatile.Write(ref isPrime, false);
break;
}
if (_kernelType == GpuKernelType.Pow2Mod)
{
if (exponent.PowModWithCycle(q, qCycle) == 1UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
else
{
UInt128 phi = q - 1UL;
if (phi <= ulong.MaxValue)
{
ulong phi64 = (ulong)phi;
if (phi64.PowModWithCycle(q, qCycle) == 1UL)
{
UInt128 halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - 1UL;
if (halfPow.BinaryGcd(q) == 1UL)
{
ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL;
ulong div = phi64.FastDiv64(exponent, divMul);
UInt128 divPow = div.PowModWithCycle(q, qCycle) - 1UL;
if (divPow.BinaryGcd(q) == 1UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
}
}
@@ -70,7 +72,7 @@ public class MersenneNumberOrderCpuTester(GpuKernelType kernelType)
auto.Next();
q = auto.CurrentQ();
}
return;
}
PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs
+7 -5
View File
@@ -1,3 +1,5 @@
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Cpu;
public class MersenneNumberResidueCpuTester
@@ -5,7 +7,7 @@ public class MersenneNumberResidueCpuTester
private ModResidueTracker? _mersenneResidueTracker;
// CPU residue variant using tracker + unrolled residue updates (no method calls in the loop).
public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
public void Scan(PrimeOrderCalculatorAccelerator gpu, ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
{
// Initialize/update Mersenne residue tracker and start a merge walk for ascending q divisors
@@ -92,7 +94,7 @@ public class MersenneNumberResidueCpuTester
// (and below) so every qualifying q reuses precomputed lengths instead of recomputing via
// tracker.MergeOrAppend.
// cycle-based quick check for small q
if (localIsPrime && tracker.MergeOrAppend(exponent, q, out isDivider) && isDivider &&
if (localIsPrime && tracker.MergeOrAppend(gpu, exponent, q, out isDivider) && isDivider &&
(q <= ulong.MaxValue ? ((ulong)q).IsPrimeCandidate() : q.IsPrimeCandidate()))
{
localIsPrime = false;
@@ -105,7 +107,7 @@ public class MersenneNumberResidueCpuTester
((allowMask >> (int)r10_1) & 1) != 0 && (r8_1 == 1UL || r8_1 == 7UL) && r3_1 != 0UL && r5_1 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
if (localIsPrime && tracker.MergeOrAppend(gpu, exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
@@ -118,7 +120,7 @@ public class MersenneNumberResidueCpuTester
((allowMask >> (int)r10_2) & 1) != 0 && (r8_2 == 1UL || r8_2 == 7UL) && r3_2 != 0UL && r5_2 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
if (localIsPrime && tracker.MergeOrAppend(gpu, exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
@@ -131,7 +133,7 @@ public class MersenneNumberResidueCpuTester
((allowMask >> (int)r10_3) & 1) != 0 && (r8_3 == 1UL || r8_3 == 7UL) && r3_3 != 0UL && r5_3 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
if (localIsPrime && tracker.MergeOrAppend(gpu, exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
PerfectNumbers.Core/DivisorCycleCache.cs
+7 -7
View File
@@ -57,13 +57,13 @@ public sealed class DivisorCycleCache
_snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
}
public ulong GetCycleLength(ulong divisor, bool skipPrimeOrderHeuristic = false)
public ulong GetCycleLength(PrimeOrderCalculatorAccelerator gpu, ulong divisor, bool skipPrimeOrderHeuristic = false)
{
if (!skipPrimeOrderHeuristic)
{
Span<ulong> result = stackalloc ulong[1];
ReadOnlySpan<ulong> singleDivisor = stackalloc ulong[1] { divisor };
GetCycleLengths(singleDivisor, result);
GetCycleLengths(gpu, singleDivisor, result);
return result[0];
}
@@ -90,12 +90,12 @@ public sealed class DivisorCycleCache
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
cycleLength = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: true);
cycleLength = MersenneDivisorCycles.CalculateCycleLength(gpu, divisor, divisorData, skipPrimeOrderHeuristic: true);
divisorPool.Return(divisorData);
return cycleLength;
}
public void GetCycleLengths(ReadOnlySpan<ulong> divisors, Span<ulong> cycles)
public void GetCycleLengths(PrimeOrderCalculatorAccelerator gpu,ReadOnlySpan<ulong> divisors, Span<ulong> cycles)
{
ulong[] snapshot = _snapshot;
int length = divisors.Length;
@@ -139,7 +139,7 @@ public sealed class DivisorCycleCache
}
else
{
ComputeCyclesCpu(divisors, cycles, missingIndices);
ComputeCyclesCpu(gpu, divisors, cycles, missingIndices);
}
}
@@ -187,7 +187,7 @@ public sealed class DivisorCycleCache
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ComputeCyclesCpu(ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices)
private static void ComputeCyclesCpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices)
{
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
for (int i = 0; i < indices.Length; i++)
@@ -195,7 +195,7 @@ public sealed class DivisorCycleCache
int targetIndex = indices[i];
ulong divisor = divisors[targetIndex];
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData);
cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(gpu, divisor, divisorData);
divisorPool.Return(divisorData);
}
}
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+38 -53
View File
@@ -20,7 +20,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
private ulong _divisorLimit;
// private bool _isConfigured;
private readonly ConcurrentBag<MersenneNumberDivisorByDivisorAccelerator>[] _resourcePool = [..AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentBag<MersenneNumberDivisorByDivisorAccelerator>())];
private readonly ConcurrentBag<MersenneNumberDivisorByDivisorAccelerator>[] _resourcePool = [.. AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentBag<MersenneNumberDivisorByDivisorAccelerator>())];
private readonly ConcurrentBag<DivisorScanSession> _sessionPool = [];
@@ -50,9 +50,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
// _isConfigured = true;
}
private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators;
public bool IsPrime(ulong prime, out bool divisorsExhausted, out ulong divisor)
public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor)
{
ulong allowedMax;
int batchCapacity;
@@ -87,7 +85,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
// GpuPrimeWorkLimiter();
// int acceleratorIndex = AcceleratorPool.Shared.Rent();
// var accelerator = _accelerators[acceleratorIndex];
var gpu = PrimeOrderCalculatorAccelerator.Rent(batchCapacity);
var acceleratorIndex = gpu.AcceleratorIndex;
// Monitor.Enter(gpuLease.ExecutionLock);
@@ -117,7 +114,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
out processedCount);
ReturnBatchResources(acceleratorIndex, resources);
PrimeOrderCalculatorAccelerator.Return(gpu);
// Monitor.Exit(gpuLease.ExecutionLock);
// GpuPrimeWorkLimiter.Release();
@@ -364,7 +360,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
{
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(gpu, divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
{
return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
return MersenneDivisorCycles.CalculateCycleLength(gpu, divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
}
return computedCycle;
@@ -481,8 +477,8 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{
private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
private readonly PrimeOrderCalculatorAccelerator _primeTesterAccelerator;
private readonly Accelerator _accelerator;
private readonly PrimeOrderCalculatorAccelerator _primeOrderCalculatorAccelerator;
private readonly AcceleratorStream _stream;
private Action<AcceleratorStream, Index1D, ArrayView<GpuDivisorPartialData>, ArrayView<int>, ArrayView<int>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>, ArrayView<int>> _kernel = null!;
private MemoryBuffer1D<GpuDivisorPartialData, Stride1D.Dense> _divisorBuffer = null!;
@@ -496,15 +492,36 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
private int _capacity;
internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner)
internal DivisorScanSession(PrimeOrderCalculatorAccelerator gpu, MersenneNumberDivisorByDivisorGpuTester owner)
{
_owner = owner;
// GpuPrimeWorkLimiter.Acquire();
_primeTesterAccelerator = PrimeOrderCalculatorAccelerator.Rent(1);
Accelerator accelerator = _primeTesterAccelerator.Accelerator;
_primeOrderCalculatorAccelerator = gpu;
Accelerator accelerator = gpu.Accelerator;
_accelerator = accelerator;
_stream = accelerator.CreateStream();
_capacity = Math.Max(1, owner._gpuBatchSize);
// lock (accelerator)
{
_divisorBuffer = accelerator.Allocate1D<GpuDivisorPartialData>(1);
_offsetBuffer = accelerator.Allocate1D<int>(1);
_countBuffer = accelerator.Allocate1D<int>(1);
_cycleBuffer = accelerator.Allocate1D<ulong>(1);
_firstHitBuffer = accelerator.Allocate1D<int>(1);
}
int desiredCapacity = PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount;
_capacity = desiredCapacity;
// lock (accelerator)
{
_exponentsBuffer = accelerator.Allocate1D<ulong>(desiredCapacity);
_hitBuffer = accelerator.Allocate1D<byte>(desiredCapacity);
}
_hostBuffer = ThreadStaticPools.UlongPool.Rent(desiredCapacity);
_kernel = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<GpuDivisorPartialData>, ArrayView<int>, ArrayView<int>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>, ArrayView<int>>(DivisorByDivisorKernels.CheckKernel)).CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView<GpuDivisorPartialData>, ArrayView<int>, ArrayView<int>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>, ArrayView<int>>>();
}
internal void Reset()
@@ -513,45 +530,16 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
private void EnsureExecutionResourcesLocked(int requiredCapacity)
{
Accelerator accelerator = _accelerator;
if (_kernel == null)
if (requiredCapacity > _capacity)
{
_kernel = KernelUtil.GetKernel(_accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<GpuDivisorPartialData>, ArrayView<int>, ArrayView<int>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>, ArrayView<int>>(DivisorByDivisorKernels.CheckKernel)).CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView<GpuDivisorPartialData>, ArrayView<int>, ArrayView<int>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>, ArrayView<int>>>();
}
_exponentsBuffer.Dispose();
_hitBuffer.Dispose();
ThreadStaticPools.UlongPool.Return(_hostBuffer, clearArray: false);
if (_divisorBuffer == null)
{
// lock (accelerator)
{
_divisorBuffer = accelerator.Allocate1D<GpuDivisorPartialData>(1);
_offsetBuffer = accelerator.Allocate1D<int>(1);
_countBuffer = accelerator.Allocate1D<int>(1);
_cycleBuffer = accelerator.Allocate1D<ulong>(1);
_firstHitBuffer = accelerator.Allocate1D<int>(1);
}
}
bool allocateBuffers = _exponentsBuffer == null;
if (!allocateBuffers && _capacity < requiredCapacity)
{
_exponentsBuffer?.Dispose();
_hitBuffer?.Dispose();
if (_hostBuffer != null)
{
ThreadStaticPools.UlongPool.Return(_hostBuffer, clearArray: false);
_hostBuffer = null!;
}
_exponentsBuffer = null!;
_hitBuffer = null!;
allocateBuffers = true;
}
if (allocateBuffers)
{
Accelerator accelerator = _accelerator;
int desiredCapacity = requiredCapacity > _capacity ? requiredCapacity : _capacity;
_capacity = desiredCapacity;
// lock (accelerator)
{
_exponentsBuffer = accelerator.Allocate1D<ulong>(desiredCapacity);
@@ -581,13 +569,12 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
ulong cycle = divisorCycle;
ulong firstPrime = primes[0];
var gpu = _primeOrderCalculatorAccelerator;
if (cycle == 0UL)
{
var gpu = _primeTesterAccelerator;
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(gpu, divisor, firstPrime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
{
cycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
cycle = MersenneDivisorCycles.CalculateCycleLength(gpu, divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
}
else
{
@@ -714,7 +701,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
return ComputeAllowedMaxDivisorGpu(prime, _divisorLimit);
}
public DivisorScanSession CreateDivisorSession()
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu)
{
// The tester is configured once at startup and the session pool relies on thread-safe collections, so the previous synchronization guard stays commented out.
// lock (_sync)
@@ -739,8 +726,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
return session;
}
return new DivisorScanSession(this);
return new DivisorScanSession(gpu, this);
}
IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() => CreateDivisorSession();
}
PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs
+3 -2
View File
@@ -1,5 +1,6 @@
using ILGPU;
using ILGPU.Runtime;
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Gpu;
@@ -86,7 +87,7 @@ public sealed class MersenneNumberDivisorGpuTester
private static (ulong divisor, uint cycle)[]? _divisorCandidates = [];
public bool IsPrime(ulong p, UInt128 d, ulong divisorCyclesSearchLimit, out bool divisorsExhausted)
public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong p, UInt128 d, ulong divisorCyclesSearchLimit, out bool divisorsExhausted)
{
ReadOnlyGpuUInt128 readOnlyDivisor;
@@ -142,7 +143,7 @@ public sealed class MersenneNumberDivisorGpuTester
continue;
}
UInt128 cycle128 = MersenneDivisorCycles.GetCycle(d);
UInt128 cycle128 = MersenneDivisorCycles.GetCycle(gpu, d);
if ((UInt128)p % cycle128 != UInt128.Zero)
{
continue;
PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs
+3 -2
View File
@@ -2,6 +2,7 @@ using System.Buffers;
using System.Collections.Concurrent;
using ILGPU;
using ILGPU.Runtime;
using PerfectNumbers.Core.Gpu.Accelerators;
using UInt128 = System.UInt128;
namespace PerfectNumbers.Core.Gpu;
@@ -131,7 +132,7 @@ public class MersenneNumberLucasLehmerGpuTester
// Configurable LL slice size to keep kernels short. Default 32.
public int SliceSize = 32;
public bool IsPrime(ulong exponent, bool runOnGpu)
public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong exponent, bool runOnGpu)
{
// Early rejections aligned with incremental/order sieves, but safe for small p:
// - If 3 | p and p != 3, then 7 | M_p -> composite.
@@ -142,7 +143,7 @@ public class MersenneNumberLucasLehmerGpuTester
return false;
}
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne(gpu))
{
return false;
}
PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs
+1 -33
View File
@@ -8,12 +8,11 @@ namespace PerfectNumbers.Core.Gpu;
public static class Pow2MontgomeryGpuCalculator
{
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong CalculateKeep(ulong exponent, in MontgomeryDivisorData divisor)
public static ulong CalculateKeep(PrimeOrderCalculatorAccelerator gpu, ulong exponent, in MontgomeryDivisorData divisor)
{
// GpuPrimeWorkLimiter.Acquire();
ulong result = 0UL;
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
int acceleratorIndex = gpu.AcceleratorIndex;
Accelerator? accelerator = gpu.Accelerator;
MemoryBuffer1D<ulong, Stride1D.Dense> exponentBuffer = gpu.Input;
@@ -32,37 +31,6 @@ public static class Pow2MontgomeryGpuCalculator
AcceleratorStreamPool.Return(acceleratorIndex, stream);
// GpuPrimeWorkLimiter.Release();
PrimeOrderCalculatorAccelerator.Return(gpu);
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong CalculateConvert(ulong exponent, in MontgomeryDivisorData divisor)
{
// GpuPrimeWorkLimiter.Acquire();
ulong result = 0UL;
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
MemoryBuffer1D<ulong, Stride1D.Dense> exponentBuffer = gpu.Input;
MemoryBuffer1D<ulong, Stride1D.Dense> resultBuffer = gpu.OutputUlong;
int acceleratorIndex = gpu.AcceleratorIndex;
AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
exponentBuffer.View.CopyFromCPU(stream, ref exponent, 1);
// We don't need to worry about any left-overs here.
// resultBuffer.MemSetToZero(stream);
var kernelLauncher = gpu.ConvertToStandardKernelLauncher!;
kernelLauncher(stream, 1, exponentBuffer.View, divisor.Modulus, divisor.NPrime, divisor.MontgomeryOne, divisor.MontgomeryTwo, divisor.MontgomeryTwoSquared, resultBuffer.View);
resultBuffer.View.CopyToCPU(stream, ref result, 1);
stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
// GpuPrimeWorkLimiter.Release();
PrimeOrderCalculatorAccelerator.Return(gpu);
return result;
}
PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
+15 -13
View File
@@ -243,7 +243,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
}
public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData divisorData)
public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData? divisorData)
{
Span<ulong> exponents = stackalloc ulong[1];
Span<ulong> remainders = stackalloc ulong[1];
@@ -254,7 +254,7 @@ internal static partial class PrimeOrderGpuHeuristics
return status;
}
public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan<ulong> exponents, ulong prime, Span<ulong> remainders, in MontgomeryDivisorData divisorData)
public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan<ulong> exponents, ulong prime, Span<ulong> remainders, in MontgomeryDivisorData? divisorData)
{
// This will never occur in production code
// if (exponents.Length == 0)
@@ -309,20 +309,20 @@ internal static partial class PrimeOrderGpuHeuristics
return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable;
}
public static GpuPow2ModStatus TryPow2Mod(in UInt128 exponent, in UInt128 prime, out UInt128 remainder)
public static GpuPow2ModStatus TryPow2Mod(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 prime, out UInt128 remainder)
{
Span<UInt128> exponents = stackalloc UInt128[1];
Span<UInt128> remainders = stackalloc UInt128[1];
exponents[0] = exponent;
GpuPow2ModStatus status = TryPow2ModBatch(exponents, prime, remainders);
GpuPow2ModStatus status = TryPow2ModBatch(gpu, exponents, prime, remainders);
remainder = remainders[0];
return status;
}
public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> remainders)
public static GpuPow2ModStatus TryPow2ModBatch(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> remainders)
{
return TryPow2ModBatchInternal(exponents, prime, remainders);
return TryPow2ModBatchInternal(gpu, exponents, prime, remainders);
}
internal static bool TryCalculateOrder(
@@ -471,7 +471,7 @@ internal static partial class PrimeOrderGpuHeuristics
return kernel;
}
private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> remainders)
private static GpuPow2ModStatus TryPow2ModBatchInternal(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> remainders)
{
if (exponents.Length == 0)
{
@@ -509,11 +509,11 @@ internal static partial class PrimeOrderGpuHeuristics
}
}
ComputeOnGpuWide(exponents, prime, target);
ComputeOnGpuWide(gpu, exponents, prime, target);
return GpuPow2ModStatus.Success;
}
private static bool TryComputeOnGpu(ReadOnlySpan<ulong> exponents, ulong prime, in MontgomeryDivisorData divisorData, Span<ulong> results)
private static bool TryComputeOnGpu(ReadOnlySpan<ulong> exponents, ulong prime, MontgomeryDivisorData? divisorData, Span<ulong> results)
{
// GpuPrimeWorkLimiter.Acquire();
int acceleratorIndex = AcceleratorPool.Shared.Rent();
@@ -540,6 +540,11 @@ internal static partial class PrimeOrderGpuHeuristics
try
{
if (divisorData is null)
{
divisorData = MontgomeryDivisorData.Empty;
}
kernelLauncher(stream, exponents.Length, exponentBuffer.View, divisorData.Modulus, divisorData.NPrime, divisorData.MontgomeryOne, divisorData.MontgomeryTwo, divisorData.MontgomeryTwoSquared, remainderBuffer.View);
remainderBuffer.View.CopyToCPU(stream, ref MemoryMarshal.GetReference(results), exponents.Length);
@@ -573,14 +578,13 @@ internal static partial class PrimeOrderGpuHeuristics
}
}
private static void ComputeOnGpuWide(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> results)
private static void ComputeOnGpuWide(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> results)
{
int length = exponents.Length;
GpuUInt128[]? rentedExponents = null;
GpuUInt128[]? rentedResults = null;
ArrayPool<GpuUInt128> gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool;
PrimeOrderCalculatorAccelerator? gpu = null;
bool poolingRequired = length > WideStackThreshold;
Span<GpuUInt128> exponentSpan = poolingRequired
? stackalloc GpuUInt128[length]
@@ -596,7 +600,6 @@ internal static partial class PrimeOrderGpuHeuristics
}
GpuUInt128 modulus = (GpuUInt128)prime;
gpu = PrimeOrderCalculatorAccelerator.Rent(1);
int acceleratorIndex = gpu.AcceleratorIndex;
Accelerator accelerator = gpu.Accelerator;
@@ -625,7 +628,6 @@ internal static partial class PrimeOrderGpuHeuristics
remainderBuffer.Dispose();
exponentBuffer.Dispose();
PrimeOrderCalculatorAccelerator.Return(gpu);
if (rentedExponents is not null)
{
PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs
+10 -10
View File
@@ -168,15 +168,15 @@ public sealed class HeuristicCombinedPrimeTester
return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10);
}
public static bool IsPrimeGpu(ulong n)
{
byte nMod10 = (byte)n.Mod10();
ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
bool isPrime = IsPrimeGpu(gpu, n, maxDivisorSquare, nMod10);
PrimeOrderCalculatorAccelerator.Return(gpu);
return isPrime;
}
// public static bool IsPrimeGpu(ulong n)
// {
// byte nMod10 = (byte)n.Mod10();
// ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
// var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
// bool isPrime = IsPrimeGpu(gpu, n, maxDivisorSquare, nMod10);
// PrimeOrderCalculatorAccelerator.Return(gpu);
// return isPrime;
// }
public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n)
{
@@ -520,7 +520,7 @@ public sealed class HeuristicCombinedPrimeTester
}
primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL;
ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(gpu, divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
cycleFromHint = false;
cycleComputed = resolvedCycle != 0UL;
return resolvedCycle;
PerfectNumbers.Core/HeuristicPrimeTester.cs
+7 -9
View File
@@ -83,17 +83,17 @@ public sealed class HeuristicPrimeTester
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeGpu(ulong n)
public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n)
{
byte nMod10 = (byte)n.Mod10();
ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
return IsPrimeGpu(n, maxDivisorSquare, nMod10);
return IsPrimeGpu(gpu, n, maxDivisorSquare, nMod10);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10) => HeuristicIsPrimeGpuCore(n, maxDivisorSquare, nMod10);
public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n, ulong maxDivisorSquare, byte nMod10) => HeuristicIsPrimeGpuCore(gpu, n, maxDivisorSquare, nMod10);
private static bool HeuristicIsPrimeGpuCore(ulong n, ulong maxDivisorSquare, byte nMod10)
private static bool HeuristicIsPrimeGpuCore(PrimeOrderCalculatorAccelerator gpu, ulong n, ulong maxDivisorSquare, byte nMod10)
{
// TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
if (maxDivisorSquare < 9UL)
@@ -106,7 +106,7 @@ public sealed class HeuristicPrimeTester
// return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB: false);
// }
bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10);
bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(gpu, n, maxDivisorSquare, nMod10);
return !compositeDetected;
}
@@ -239,10 +239,9 @@ public sealed class HeuristicPrimeTester
}
private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
private static bool HeuristicTrialDivisionGpuDetectsDivisor(PrimeOrderCalculatorAccelerator gpu, ulong n, ulong maxDivisorSquare, byte nMod10)
{
// GpuPrimeWorkLimiter.Acquire();
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
int acceleratorIndex = gpu.AcceleratorIndex;
var stream = AcceleratorStreamPool.Rent(acceleratorIndex);
var flagView1D = gpu.OutputInt!.View;
@@ -313,7 +312,6 @@ public sealed class HeuristicPrimeTester
// }
AcceleratorStreamPool.Return(acceleratorIndex, stream);
PrimeOrderCalculatorAccelerator.Return(gpu);
// GpuPrimeWorkLimiter.Release();
return compositeDetected;
}
@@ -487,7 +485,7 @@ public sealed class HeuristicPrimeTester
}
primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL;
ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(gpu, divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
cycleFromHint = false;
cycleComputed = resolvedCycle != 0UL;
return resolvedCycle;
PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs
+4 -2
View File
@@ -1,3 +1,5 @@
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core;
public interface IMersenneNumberDivisorByDivisorTester
@@ -10,9 +12,9 @@ public interface IMersenneNumberDivisorByDivisorTester
ulong GetAllowedMaxDivisor(ulong prime);
bool IsPrime(ulong prime, out bool divisorsExhausted, out ulong divisor);
bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor);
IDivisorScanSession CreateDivisorSession();
IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu);
void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues);
PerfectNumbers.Core/MersenneDivisorCycles.cs
+14 -9
View File
@@ -75,7 +75,7 @@ public class MersenneDivisorCycles
}
}
if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL)
if (TryCalculateCycleLengthHeuristic(gpu, divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL)
{
return heuristicCycle == exponent;
}
@@ -163,7 +163,7 @@ public class MersenneDivisorCycles
}
}
public ulong GetCycle(ulong divisor, in MontgomeryDivisorData divisorData)
public ulong GetCycle(PrimeOrderCalculatorAccelerator gpu, ulong divisor, in MontgomeryDivisorData divisorData)
{
// Fast-path: in-memory array for small divisors
if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
@@ -208,10 +208,10 @@ public class MersenneDivisorCycles
// TODO: Route this miss to an ephemeral single-cycle computation on the device selected by the current
// configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot
// so we respect the large-p memory limits and avoid extra cache locks.
return CalculateCycleLength(divisor, divisorData);
return CalculateCycleLength(gpu, divisor, divisorData);
}
public static UInt128 GetCycle(UInt128 divisor)
public static UInt128 GetCycle(PrimeOrderCalculatorAccelerator gpu, UInt128 divisor)
{
// For divisor = 2^k, cycle is 1
if ((divisor & (divisor - UInt128.One)) == UInt128.Zero)
@@ -224,12 +224,13 @@ public class MersenneDivisorCycles
ulong divisor64 = (ulong)divisor;
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor64);
ulong cycle = CalculateCycleLength(divisor64, divisorData);
ulong cycle = CalculateCycleLength(gpu, divisor64, divisorData);
divisorPool.Return(divisorData);
return (UInt128)cycle;
}
UInt128 wideOrder = PrimeOrderCalculator.Calculate(
gpu,
divisor,
previousOrder: null,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
@@ -285,6 +286,7 @@ public class MersenneDivisorCycles
// }
ulong computedOrder = PrimeOrderCalculator.Calculate(
gpu,
divisor,
previousOrder: null,
divisorData,
@@ -743,6 +745,7 @@ public class MersenneDivisorCycles
// (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength));
int localCycleIndex = 0;
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
for (ulong divisor = localStart; divisor <= localEnd; divisor++)
{
if ((divisor & 1UL) == 0UL)
@@ -762,11 +765,12 @@ public class MersenneDivisorCycles
// unrolled-hex calculator so the CPU generator matches
// the GPU benchmark leader for large divisors.
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
ulong cycle = CalculateCycleLength(divisor, divisorData);
ulong cycle = CalculateCycleLength(gpu, divisor, divisorData);
divisorPool.Return(divisorData);
localCycles[localCycleIndex++] = (divisor, cycle);
}
PrimeOrderCalculatorAccelerator.Return(gpu);
if (localCycleIndex == 0)
{
ThreadStaticPools.Return(localCycles);
@@ -929,9 +933,9 @@ public class MersenneDivisorCycles
public static ulong CalculateCycleLengthGpu(ulong divisor) => DivisorCycleKernels.CalculateCycleLengthGpu(divisor);
public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false)
public static ulong CalculateCycleLength(PrimeOrderCalculatorAccelerator gpu, ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false)
{
if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong cycleLength, skipPrimeOrderHeuristic))
if (TryCalculateCycleLengthHeuristic(gpu, divisor, divisorData, out ulong cycleLength, skipPrimeOrderHeuristic))
return cycleLength;
return CalculateCycleLengthFallback(divisor);
@@ -940,7 +944,7 @@ public class MersenneDivisorCycles
// [ThreadStatic]
// private static HeuristicPrimeTester? _primeTester;
internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false)
internal static bool TryCalculateCycleLengthHeuristic(PrimeOrderCalculatorAccelerator gpu, ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false)
{
if ((divisor & (divisor - 1UL)) == 0UL)
{
@@ -966,6 +970,7 @@ public class MersenneDivisorCycles
// if (!skipPrimeOrderHeuristic && HeuristicCombinedPrimeTester.IsPrimeGpu(divisor))
{
ulong computedOrder = PrimeOrderCalculator.Calculate(
gpu,
divisor,
previousOrder: null,
divisorData,
PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs
+9 -4
View File
@@ -239,10 +239,10 @@ public static class MersenneNumberDivisorByDivisorTester
int chunkSize = primesPerTask <= 0 ? 1 : primesPerTask;
void ProcessPrime(ulong prime)
void ProcessPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime)
{
Console.WriteLine($"Processing {prime}");
bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted, out ulong divisor);
bool isPrime = tester.IsPrime(gpu, prime, out bool divisorsExhausted, out ulong divisor);
if (!isPrime)
{
@@ -259,10 +259,13 @@ public static class MersenneNumberDivisorByDivisorTester
if (workerCount == 1)
{
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
foreach (ulong prime in filteredPrimes)
{
ProcessPrime(prime);
ProcessPrime(gpu, prime);
}
PrimeOrderCalculatorAccelerator.Return(gpu);
}
else
{
@@ -295,11 +298,13 @@ public static class MersenneNumberDivisorByDivisorTester
startGate.Wait();
Console.WriteLine($"Task started for range {rangeStart}");
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
for (int index = rangeStart; index < rangeEnd; index++)
{
ProcessPrime(filteredPrimes[index]);
ProcessPrime(gpu, filteredPrimes[index]);
}
PrimeOrderCalculatorAccelerator.Return(gpu);
Console.WriteLine($"Task finished for range {rangeStart}");
PrimeOrderCalculatorAccelerator.DisposeAll();
},
PerfectNumbers.Core/MersenneNumberTester.cs
+13 -10
View File
@@ -63,6 +63,7 @@ public sealed class MersenneNumberTester(
UInt128 twoP = (UInt128)exponent << 1; // 2 * p
LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
// CPU path: compute orders directly and cache
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
if (!_useGpuOrder)
{
UInt128 k = UInt128.One;
@@ -97,10 +98,11 @@ public sealed class MersenneNumberTester(
continue;
}
ulong ord = q.CalculateOrder();
ulong ord = q.CalculateOrder(gpu);
_orderCache[q] = ord;
}
PrimeOrderCalculatorAccelerator.Return(gpu);
return;
}
@@ -190,7 +192,7 @@ public sealed class MersenneNumberTester(
ulong order = orders[i];
if (order == 0UL)
{
order = qs[offset + i].CalculateOrder();
order = qs[offset + i].CalculateOrder(gpu);
}
if (useOrderCache)
@@ -205,11 +207,12 @@ public sealed class MersenneNumberTester(
}
AcceleratorStreamPool.Return(acceleratorIndex,stream);
PrimeOrderCalculatorAccelerator.Return(gpu);
// GpuPrimeWorkLimiter.Release();
uInt128Pool.Return(qs);
}
public bool IsMersennePrime(ulong exponent)
public bool IsMersennePrime(PrimeOrderCalculatorAccelerator gpu, ulong exponent)
{
// Safe early rejections using known orders for tiny primes:
// 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime.
@@ -220,7 +223,7 @@ public sealed class MersenneNumberTester(
return false;
}
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne(gpu))
{
return false;
}
@@ -251,7 +254,7 @@ public sealed class MersenneNumberTester(
}
else
{
_residueCpuTester!.Scan(exponent, twoP, lastDigit, maxK, ref prePrime);
_residueCpuTester!.Scan(gpu, exponent, twoP, lastDigit, maxK, ref prePrime);
}
return prePrime;
@@ -270,7 +273,7 @@ public sealed class MersenneNumberTester(
}
else
{
_incrementalCpuTester.Scan(exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime);
_incrementalCpuTester.Scan(gpu, exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime);
}
if (!prePrime)
@@ -279,8 +282,8 @@ public sealed class MersenneNumberTester(
}
return _useGpuLucas
? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
: _lucasLehmerCpuTester.IsPrime(exponent);
? _lucasLehmerGpuTester.IsPrime(gpu, exponent, runOnGpu: true)
: _lucasLehmerCpuTester.IsPrime(gpu, exponent);
}
bool isPrime = true;
@@ -293,7 +296,7 @@ public sealed class MersenneNumberTester(
}
else
{
_orderCpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime);
_orderCpuTester.Scan(gpu, exponent, twoP, lastDigit, maxK, ref isPrime);
}
}
else
@@ -304,7 +307,7 @@ public sealed class MersenneNumberTester(
}
else
{
_incrementalCpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime);
_incrementalCpuTester.Scan(gpu, exponent, twoP, lastDigit, maxK, ref isPrime);
}
}
PerfectNumbers.Core/ModResidueTracker.cs
+8 -8
View File
@@ -1,5 +1,5 @@
using System.Runtime.CompilerServices;
using PerfectNumbers.Core.Gpu.Accelerators;
using UInt128 = System.UInt128;
namespace PerfectNumbers.Core;
@@ -39,7 +39,7 @@ public sealed class ModResidueTracker
/// (number must be non-decreasing across calls). Adds and initializes new
/// divisors on first use.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public bool IsDivisible(UInt128 number, UInt128 divisor)
public bool IsDivisible(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor)
{
if (divisor == UInt128.Zero)
{
@@ -52,7 +52,7 @@ public sealed class ModResidueTracker
int idx = LowerBound(divisors, divisor);
if (idx == divisors.Count || divisors[idx] != divisor)
{
UInt128 residue = ComputeInitialResidue(number, divisor);
UInt128 residue = ComputeInitialResidue(gpu, number, divisor);
divisors.Insert(idx, divisor);
_residues.Insert(idx, residue);
if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended)
@@ -81,7 +81,7 @@ public sealed class ModResidueTracker
// Merge or append current candidate divisor (ascending order required across calls).
// Returns true and sets 'divisible' based on residue at 'number'.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public bool MergeOrAppend(UInt128 number, UInt128 divisor, out bool divisible)
public bool MergeOrAppend(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor, out bool divisible)
{
// Assume number is non-decreasing; advance residues lazily when number changes.
List<UInt128> divisors = _divisors;
@@ -97,7 +97,7 @@ public sealed class ModResidueTracker
List<UInt128> residues = _residues;
if (divisors.Count == 0 || (_hasLastAppended && divisor >= _lastAppended && divisors[^1] <= divisor))
{
residue = ComputeInitialResidue(number, divisor);
residue = ComputeInitialResidue(gpu, number, divisor);
divisors.Add(divisor);
residues.Add(residue);
_hasLastAppended = true;
@@ -123,7 +123,7 @@ public sealed class ModResidueTracker
// Insert at cursor to keep order
residue = ComputeInitialResidue(number, divisor);
residue = ComputeInitialResidue(gpu, number, divisor);
divisors.Insert(cursor, divisor);
residues.Insert(cursor, residue);
divisible = residue == UInt128.Zero;
@@ -238,7 +238,7 @@ public sealed class ModResidueTracker
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private UInt128 ComputeInitialResidue(UInt128 number, UInt128 divisor)
private UInt128 ComputeInitialResidue(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor)
{
if (_model == ResidueModel.Identity)
{
@@ -266,7 +266,7 @@ public sealed class ModResidueTracker
if (pow != UInt128.One)
return (pow - UInt128.One) % divisor;
ulong ord = divisor.CalculateOrder();
ulong ord = divisor.CalculateOrder(gpu);
if (ord == 0UL || number % (UInt128)ord != UInt128.Zero)
return UInt128.One; // does NOT divide
PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
+1 -5
View File
@@ -44,7 +44,7 @@ internal static partial class PrimeOrderCalculator
}
}
private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
private static ulong CalculateInternal(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
{
// TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31.
if (prime <= 3UL)
@@ -54,10 +54,8 @@ internal static partial class PrimeOrderCalculator
ulong phi = prime - 1UL;
PrimeOrderCalculatorAccelerator gpu = PrimeOrderCalculatorAccelerator.Rent(1);
if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(gpu, prime, previousOrder, config, divisorData, out ulong gpuOrder))
{
PrimeOrderCalculatorAccelerator.Return(gpu);
return gpuOrder;
}
@@ -73,13 +71,11 @@ internal static partial class PrimeOrderCalculator
// result = CalculateByFactorizationCpu(gpu, prime, divisorData);
phiFactors.Dispose();
PrimeOrderCalculatorAccelerator.Return(gpu);
return result;
}
result = RunHeuristicPipelineCpu(gpu, prime, previousOrder, config, divisorData, phi, phiFactors);
phiFactors.Dispose();
PrimeOrderCalculatorAccelerator.Return(gpu);
return result;
}
PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs
+32 -29
View File
@@ -2,6 +2,7 @@ using System.Buffers;
using System.Diagnostics;
using System.Numerics;
using PerfectNumbers.Core.Gpu;
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core;
@@ -82,7 +83,7 @@ internal static partial class PrimeOrderCalculator
Console.WriteLine(messageFactory());
}
}
private static UInt128 CalculateWideInternal(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData? divisorData, in PrimeOrderSearchConfig config)
private static UInt128 CalculateWideInternal(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData? divisorData, in PrimeOrderSearchConfig config)
{
if (prime <= UInt128.One)
{
@@ -100,11 +101,11 @@ internal static partial class PrimeOrderCalculator
UInt128 result;
if (phiFactors.Factors is null)
{
result = FinishStrictlyWide(prime, divisorData);
result = FinishStrictlyWide(gpu, prime, divisorData);
}
else
{
result = RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
result = RunHeuristicPipelineWide(gpu, prime, previousOrder, divisorData, config, phi, phiFactors);
}
debugScope.Dispose();
@@ -112,32 +113,33 @@ internal static partial class PrimeOrderCalculator
}
private static UInt128 RunHeuristicPipelineWide(
PrimeOrderCalculatorAccelerator gpu,
in UInt128 prime,
in UInt128? previousOrder,
in MontgomeryDivisorData divisorData,
in MontgomeryDivisorData? divisorData,
in PrimeOrderSearchConfig config,
in UInt128 phi,
in PartialFactorResult128 phiFactors)
{
if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
if (phiFactors.FullyFactored && TrySpecialMaxWide(gpu, phi, prime, divisorData, phiFactors))
{
return phi;
}
UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
candidateOrder = ExponentLoweringWide(candidateOrder, prime, divisorData, phiFactors);
UInt128 candidateOrder = InitializeStartingOrderWide(gpu, prime, phi, divisorData);
candidateOrder = ExponentLoweringWide(gpu, candidateOrder, prime, divisorData, phiFactors);
if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
if (TryConfirmOrderWide(gpu, prime, candidateOrder, divisorData, config))
{
return candidateOrder;
}
if (config.Mode == PrimeOrderMode.Strict)
{
return FinishStrictlyWide(prime, divisorData);
return FinishStrictlyWide(gpu, prime, divisorData);
}
if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
if (TryHeuristicFinishWide(gpu, prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
{
return order;
}
@@ -145,12 +147,12 @@ internal static partial class PrimeOrderCalculator
return candidateOrder;
}
private static UInt128 FinishStrictlyWide(in UInt128 prime, in MontgomeryDivisorData divisorData)
private static UInt128 FinishStrictlyWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData? divisorData)
{
return CalculateByFactorizationWide(prime, divisorData);
return CalculateByFactorizationWide(gpu, prime, divisorData);
}
private static bool TrySpecialMaxWide(in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
private static bool TrySpecialMaxWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors)
{
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
int length = factors.Count;
@@ -158,7 +160,7 @@ internal static partial class PrimeOrderCalculator
{
UInt128 factor = factorSpan[i].Value;
UInt128 reduced = phi / factor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One)
{
return false;
}
@@ -167,14 +169,14 @@ internal static partial class PrimeOrderCalculator
return true;
}
private static UInt128 InitializeStartingOrderWide(in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData divisorData)
private static UInt128 InitializeStartingOrderWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData? divisorData)
{
UInt128 order = phi;
UInt128 mod8 = prime & (UInt128)7UL;
if (mod8 == UInt128.One || mod8 == (UInt128)7UL)
{
UInt128 half = phi >> 1;
if (Pow2ModWide(half, prime, divisorData) == UInt128.One)
if (Pow2ModWide(gpu, half, prime, divisorData) == UInt128.One)
{
order = half;
}
@@ -183,7 +185,7 @@ internal static partial class PrimeOrderCalculator
return order;
}
private static UInt128 ExponentLoweringWide(UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
private static UInt128 ExponentLoweringWide(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
{
ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool;
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
@@ -215,7 +217,7 @@ internal static partial class PrimeOrderCalculator
if ((order % primeFactor) == UInt128.Zero)
{
UInt128 reduced = order / primeFactor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One)
{
order = reduced;
continue;
@@ -230,14 +232,14 @@ internal static partial class PrimeOrderCalculator
return order;
}
private static bool TryConfirmOrderWide(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
private static bool TryConfirmOrderWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
{
if (order == UInt128.Zero)
{
return false;
}
if (Pow2ModWide(order, prime, divisorData) != UInt128.One)
if (Pow2ModWide(gpu, order, prime, divisorData) != UInt128.One)
{
return false;
}
@@ -277,7 +279,7 @@ internal static partial class PrimeOrderCalculator
}
reduced /= primeFactor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One)
{
return false;
}
@@ -288,6 +290,7 @@ internal static partial class PrimeOrderCalculator
}
private static bool TryHeuristicFinishWide(
PrimeOrderCalculatorAccelerator gpu,
in UInt128 prime,
in UInt128 order,
in UInt128? previousOrder,
@@ -347,12 +350,12 @@ internal static partial class PrimeOrderCalculator
UInt128 candidate = candidates[i];
powUsed++;
if (Pow2ModWide(candidate, prime, divisorData) != UInt128.One)
if (Pow2ModWide(gpu, candidate, prime, divisorData) != UInt128.One)
{
continue;
}
if (!TryConfirmCandidateWide(prime, candidate, config, ref powUsed, powBudget, divisorData))
if (!TryConfirmCandidateWide(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData))
{
continue;
}
@@ -536,7 +539,7 @@ internal static partial class PrimeOrderCalculator
}
}
private static bool TryConfirmCandidateWide(in UInt128 prime, in UInt128 candidate, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData)
private static bool TryConfirmCandidateWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData)
{
PartialFactorResult128 factorization = PartialFactorWide(candidate, config);
if (factorization.Factors is null)
@@ -579,7 +582,7 @@ internal static partial class PrimeOrderCalculator
}
powUsed++;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
if (Pow2ModWide(gpu, reduced, prime, divisorData) == UInt128.One)
{
return false;
}
@@ -855,7 +858,7 @@ internal static partial class PrimeOrderCalculator
return 64 + BitOperations.TrailingZeroCount(high);
}
private static UInt128 CalculateByFactorizationWide(in UInt128 prime, in MontgomeryDivisorData divisorData)
private static UInt128 CalculateByFactorizationWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData? divisorData)
{
UInt128 phi = prime - UInt128.One;
Dictionary<UInt128, int> counts = new(capacity: 8);
@@ -879,7 +882,7 @@ internal static partial class PrimeOrderCalculator
if ((order % primeFactor) != UInt128.Zero)
{
UInt128 candidate = order / primeFactor;
if (Pow2ModWide(candidate, prime, divisorData) == UInt128.One)
if (Pow2ModWide(gpu, candidate, prime, divisorData) == UInt128.One)
{
order = candidate;
}
@@ -1026,7 +1029,7 @@ internal static partial class PrimeOrderCalculator
return true;
}
private static UInt128 Pow2ModWide(in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData)
private static UInt128 Pow2ModWide(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData? divisorData)
{
if (modulus == UInt128.One)
{
@@ -1046,7 +1049,7 @@ internal static partial class PrimeOrderCalculator
}
}
GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, modulus, out UInt128 wideRemainder);
GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, modulus, out UInt128 wideRemainder);
if (wideStatus == GpuPow2ModStatus.Success)
{
return wideRemainder;
PerfectNumbers.Core/PrimeOrderCalculator.cs
+6 -3
View File
@@ -1,4 +1,5 @@
using System.Runtime.CompilerServices;
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core;
@@ -48,6 +49,7 @@ internal static partial class PrimeOrderCalculator
// TODO: Remove branching on the CPU / GPU paths
// TODO: Split big non-static functions into smaller / extract static code to limit JIT time
public static ulong Calculate(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
@@ -55,12 +57,13 @@ internal static partial class PrimeOrderCalculator
PrimeOrderHeuristicDevice device)
{
var scope = UsePow2Mode(device);
ulong order = CalculateInternal(prime, previousOrder, divisorData, config);
ulong order = CalculateInternal(gpu, prime, previousOrder, divisorData, config);
scope.Dispose();
return order;
}
public static UInt128 Calculate(
PrimeOrderCalculatorAccelerator gpu,
in UInt128 prime,
in UInt128? previousOrder,
in PrimeOrderSearchConfig config,
@@ -88,14 +91,14 @@ internal static partial class PrimeOrderCalculator
ulong prime64 = (ulong)prime;
Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
divisorData = divisorPool.FromModulus(prime64);
ulong order64 = Calculate(prime64, previous, divisorData, config, device);
ulong order64 = Calculate(gpu, prime64, previous, divisorData, config, device);
divisorPool.Return(divisorData);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
}
else
{
divisorData = MontgomeryDivisorData.Empty;
result = CalculateWideInternal(prime, previousOrder, divisorData, config);
result = CalculateWideInternal(gpu, prime, previousOrder, divisorData, config);
}
scope.Dispose();
PerfectNumbers.Core/PrimeTester.cs
+34 -38
View File
@@ -87,45 +87,45 @@ public sealed class PrimeTester
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeGpu(ulong n)
{
// [MethodImpl(MethodImplOptions.AggressiveInlining)]
// public static bool IsPrimeGpu(ulong n)
// {
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
int acceleratorIndex = gpu.AcceleratorIndex;
var inputView = gpu.Input.View;
var outputView = gpu.OutputByte.View;
// var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
// int acceleratorIndex = gpu.AcceleratorIndex;
// var inputView = gpu.Input.View;
// var outputView = gpu.OutputByte.View;
// GpuPrimeWorkLimiter.Acquire();
AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
inputView.CopyFromCPU(stream, ref n, 1);
// // GpuPrimeWorkLimiter.Acquire();
// AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
// inputView.CopyFromCPU(stream, ref n, 1);
var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher;
// var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher;
kernelLauncher(
stream,
1,
inputView,
gpu.DevicePrimesLastOne,
gpu.DevicePrimesLastSeven,
gpu.DevicePrimesLastThree,
gpu.DevicePrimesLastNine,
gpu.DevicePrimesPow2LastOne,
gpu.DevicePrimesPow2LastSeven,
gpu.DevicePrimesPow2LastThree,
gpu.DevicePrimesPow2LastNine,
outputView);
// kernelLauncher(
// stream,
// 1,
// inputView,
// gpu.DevicePrimesLastOne,
// gpu.DevicePrimesLastSeven,
// gpu.DevicePrimesLastThree,
// gpu.DevicePrimesLastNine,
// gpu.DevicePrimesPow2LastOne,
// gpu.DevicePrimesPow2LastSeven,
// gpu.DevicePrimesPow2LastThree,
// gpu.DevicePrimesPow2LastNine,
// outputView);
byte flag = 0;
outputView.CopyToCPU(stream, ref flag, 1);
stream.Synchronize();
// byte flag = 0;
// outputView.CopyToCPU(stream, ref flag, 1);
// stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
PrimeOrderCalculatorAccelerator.Return(gpu);
// GpuPrimeWorkLimiter.Release();
// AcceleratorStreamPool.Return(acceleratorIndex, stream);
// PrimeOrderCalculatorAccelerator.Return(gpu);
// // GpuPrimeWorkLimiter.Release();
return flag != 0;
}
// return flag != 0;
// }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n)
@@ -189,10 +189,9 @@ public sealed class PrimeTester
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void IsPrimeBatchGpu(ReadOnlySpan<ulong> values, Span<byte> results)
public static void IsPrimeBatchGpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> values, Span<byte> results)
{
// GpuPrimeWorkLimiter.Acquire();
var gpu = PrimeOrderCalculatorAccelerator.Rent(GpuBatchSize);
int acceleratorIndex = gpu.AcceleratorIndex;
int totalLength = values.Length;
int batchSize = GpuBatchSize;
@@ -236,7 +235,6 @@ public sealed class PrimeTester
stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
PrimeOrderCalculatorAccelerator.Return(gpu);
// GpuPrimeWorkLimiter.Release();
}
@@ -265,7 +263,7 @@ public sealed class PrimeTester
return n.BinaryGcd(m) != 1UL;
}
internal static void SharesFactorWithMaxExponentBatch(ReadOnlySpan<ulong> values, Span<byte> results)
internal static void SharesFactorWithMaxExponentBatch(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> values, Span<byte> results)
{
// TODO: Route this batch helper through the shared GPU kernel pool from
// GpuUInt128BinaryGcdBenchmarks so we reuse cached kernels, pinned host buffers,
@@ -275,7 +273,6 @@ public sealed class PrimeTester
int length = values.Length;
var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
int acceleratorIndex = gpu.AcceleratorIndex;
gpu.EnsureCapacity(0, length);
MemoryBuffer1D<ulong, Stride1D.Dense>? inputBuffer = gpu.Input;
@@ -291,7 +288,6 @@ public sealed class PrimeTester
stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
PrimeOrderCalculatorAccelerator.Return(gpu);
}
}
PerfectNumbers.Core/UInt128Extensions.cs
+732 -731
View File
File diff suppressed because it is too large. Load diff
PerfectNumbers.Core/ULongExtensions.Gpu.cs
+6 -10
View File
@@ -44,27 +44,23 @@ public static partial class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWindowedKeepGpu(this ulong exponent, in MontgomeryDivisorData divisor)
=> Pow2MontgomeryGpuCalculator.CalculateKeep(exponent, divisor);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWindowedConvertGpu(this ulong exponent, in MontgomeryDivisorData divisor)
=> Pow2MontgomeryGpuCalculator.CalculateConvert(exponent, divisor);
public static ulong Pow2MontgomeryModWindowedKeepGpu(this ulong exponent, PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisor)
=> Pow2MontgomeryGpuCalculator.CalculateKeep(gpu, exponent, divisor);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWindowedConvertGpu(this ulong exponent, PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisor)
=> Pow2MontgomeryGpuCalculator.CalculateConvert(gpu, exponent, divisor);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, PrimeOrderCalculatorAccelerator gpu, ulong cycleLength, in MontgomeryDivisorData divisor)
{
ulong rotationCount = exponent % cycleLength;
return Pow2MontgomeryGpuCalculator.CalculateConvert(rotationCount, divisor);
return Pow2MontgomeryGpuCalculator.CalculateConvert(gpu, rotationCount, divisor);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisor)
{
return Pow2MontgomeryGpuCalculator.CalculateConvert(reducedExponent, divisor);
return Pow2MontgomeryGpuCalculator.CalculateConvert(gpu, reducedExponent, divisor);
}
}
PerfectNumbers.Core/ULongExtensions.cs
+8 -7
View File
@@ -1,5 +1,6 @@
using System.Numerics;
using System.Runtime.CompilerServices;
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core;
@@ -98,7 +99,7 @@ public static partial class ULongExtensions
}
}
public static ulong CalculateOrder(this ulong q)
public static ulong CalculateOrder(this ulong q, PrimeOrderCalculatorAccelerator gpu)
{
if (q <= 2UL)
{
@@ -111,7 +112,7 @@ public static partial class ULongExtensions
int i = 0, primesLength = smallPrimes.Length;
UInt128 q128 = q,
cycle = MersenneDivisorCycles.GetCycle(q128);
cycle = MersenneDivisorCycles.GetCycle(gpu, q128);
// TODO: When the shared cycle snapshot cannot serve this divisor, trigger an on-demand
// GPU computation (respecting the configured device) without promoting the result into
// the cache so the order calculator still benefits from cycle stepping while keeping the
@@ -860,7 +861,7 @@ public static partial class ULongExtensions
/// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus.
/// Only for modulus >= 10 and reasonable size.
/// </summary>
public static UInt128 PowModCrt(this ulong exponent, UInt128 modulus, MersenneDivisorCycles cycles)
public static UInt128 PowModCrt(this ulong exponent, PrimeOrderCalculatorAccelerator gpu, UInt128 modulus, MersenneDivisorCycles cycles)
{
if (modulus < 10)
return PowMod(exponent, modulus); // fallback to classic
@@ -875,7 +876,7 @@ public static partial class ULongExtensions
for (; modulusCandidate <= modulus; modulusCandidate++)
{
cycle = MersenneDivisorCycles.GetCycle(modulusCandidate);
cycle = MersenneDivisorCycles.GetCycle(gpu, modulusCandidate);
remainderForCandidate = cycle > zero
? PowModWithCycle(exponent, modulusCandidate, cycle)
: PowMod(exponent, modulusCandidate);
@@ -935,7 +936,7 @@ public static partial class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool SharesFactorWithExponentMinusOne(this ulong exponent)
public static bool SharesFactorWithExponentMinusOne(this ulong exponent, PrimeOrderCalculatorAccelerator gpu)
{
ulong prime, value = exponent - 1UL;
value >>= BitOperations.TrailingZeroCount(value);
@@ -951,7 +952,7 @@ public static partial class ULongExtensions
continue;
}
if (exponent % prime.CalculateOrder() == 0UL)
if (exponent % prime.CalculateOrder(gpu) == 0UL)
{
return true;
}
@@ -963,7 +964,7 @@ public static partial class ULongExtensions
while (value % prime == 0UL);
}
if (value > 1UL && exponent % value.CalculateOrder() == 0UL)
if (value > 1UL && exponent % value.CalculateOrder(gpu) == 0UL)
{
return true;
}