Compare commits

...
Author SHA1 Message Date
Marek Lemańczyk 6027a73839 Split GPU heuristic tables by scope 2025-11-01 09:22:44 +01:00
Marek Lemańczyk 936ce239db Merge pull request #251 from mlemanczyk/codex/refactor-prime-checking-on-gpu-n5t2pm
Warm GPU resources before thread startup
2025-11-01 04:09:00 +01:00
mlemanczyk 22f20bfe98 Various small optimizations 2025-11-01 03:42:36 +01:00
Marek Lemańczyk b7c34f4960 Preload divisor generator tables for GPU pool contexts 2025-10-31 23:01:28 +01:00
Marek Lemańczyk 469c8057d3 Flatten heuristic trial division bounds checks 2025-10-31 22:45:59 +01:00
Marek Lemańczyk 8c662c5c85 Address feedback on GPU pooling and heuristic factoring 2025-10-31 21:47:42 +01:00
Marek Lemańczyk 567b92792c Use span growth instead of dictionaries in partial factor 2025-10-31 18:52:00 +01:00
Marek Lemańczyk 6864b316d5 Warm GPU resources before thread startup 2025-10-31 16:56:25 +01:00
mlemanczyk e14b77e504 Batch size made dynamic - dependent on the value 2025-10-30 02:43:15 +01:00
Marek Lemańczyk 30af16a4f6 Merge pull request #249 from mlemanczyk/codex/optimize-heuristiccombinedprimetester-calculations
Use divisor square tables in heuristic combined prime tester
2025-10-29 22:31:52 +01:00
mlemanczyk fc5700f82c IsPrime with cancellation token removed. 2025-10-29 22:31:28 +01:00
Marek Lemańczyk f9504a6a19 Use atomic flag writes in heuristic trial division 2025-10-29 20:40:40 +01:00
Marek Lemańczyk 5ac05489c1 Use precomputed divisor squares in heuristic prime tester 2025-10-29 20:22:41 +01:00
Marek Lemańczyk f7a0627ce5 Align heuristic prime testers with squared limits 2025-10-29 19:35:59 +01:00
Marek Lemańczyk 5536dc4865 Fix GPU by-divisor hang for small primes 2025-10-29 19:21:41 +01:00
Marek Lemańczyk c61df1c33b Inline batch copies for GPU prime tester 2025-10-29 19:21:32 +01:00
Marek Lemańczyk 2b0fc5be13 Pool reusable GPU leases with embedded kernel state 2025-10-29 18:22:40 +01:00
Marek Lemańczyk a4ac82e74d Inline GPU kernel state into prime tester leases 2025-10-29 18:13:51 +01:00
Marek Lemańczyk 23fa56b356 Embed kernel state lookup in GPU leases 2025-10-29 18:13:46 +01:00
Marek Lemańczyk 08499fdda1 Refine GPU prime tester leasing 2025-10-29 17:39:38 +01:00
Marek Lemańczyk 6395b06b19 Optimize single GPU prime path and reuse lease buffers 2025-10-29 17:20:21 +01:00
Marek Lemańczyk 83a6c6c558 Scope GPU prime tester buffers to leases 2025-10-29 17:05:17 +01:00
Marek Lemańczyk ca53ddc798 Inline GPU lease buffers into the context lease 2025-10-29 16:53:29 +01:00
Marek Lemańczyk aab785650d Manage GPU lease buffers without kernel scratch pool 2025-10-29 16:53:22 +01:00
Marek Lemańczyk 13b59abb1c Embed GPU scratch buffers in context leases 2025-10-29 16:25:38 +01:00
Marek Lemańczyk b753654ebb Enforce valid divisor selectors in heuristic combined tester 2025-10-29 15:32:40 +01:00
Marek Lemańczyk b0158a844a Use divisor squares to bound heuristic combined trial division 2025-10-29 15:13:54 +01:00
mlemanczyk 4ab77ba170 Merge branch 'codex/analyze-gc-instances-in-checkdivisors-method' of https://github.com/mlemanczyk/even-perfect-numbers-scanner into codex/analyze-gc-instances-in-checkdivisors-method 2025-10-29 14:06:08 +01:00
Marek Lemańczyk 47eee81d3c Merge pull request #248 from mlemanczyk/codex/optimize-heuristicprimetester-for-divisor-selection
Add filtered group B divisor tables for heuristic GPU checks
2025-10-29 14:05:02 +01:00
mlemanczyk b2865b54ce Tiny optimization 2025-10-29 14:04:25 +01:00
mlemanczyk 4c79c4087e Various optimizations 2025-10-29 14:00:38 +01:00
mlemanczyk 3147728d5b Benchmarks for heuristic combined methods added 2025-10-29 04:53:13 +01:00
mlemanczyk 880b1cd962 Revert "Unroll divisor-by-divisor GPU kernel loops"
This reverts commit 5d14efb739.
2025-10-29 04:51:28 +01:00
Marek Lemańczyk 5d14efb739 Unroll divisor-by-divisor GPU kernel loops 2025-10-29 04:31:45 +01:00
Marek Lemańczyk f49963c4f6 Restore Math.Sqrt heuristic bound 2025-10-29 01:43:12 +01:00
Marek Lemańczyk 5a9e103bd5 Replace heuristic sqrt with log2-based estimate 2025-10-29 01:29:12 +01:00
Marek Lemańczyk eeeff42be6 Replace heuristic sqrt with integer algorithm 2025-10-29 01:02:31 +01:00
Marek Lemańczyk b54f4b6bd6 Simplify combined heuristic divisor selection 2025-10-29 00:44:46 +01:00
Marek Lemańczyk b2cb17a96f Return to one-A-one-B combined schedule 2025-10-29 00:44:39 +01:00
Marek Lemańczyk f993561cd5 Eagerly build combined divisors and restore two-A-one-B default 2025-10-28 23:43:15 +01:00
Marek Lemańczyk e8d163f595 Add three-A-two-B combined divisor schedule 2025-10-28 23:22:16 +01:00
Marek Lemańczyk 465b40bf35 Add three-A-one-B schedule to combined tester 2025-10-28 23:09:41 +01:00
Marek Lemańczyk 4c445f2a42 Switch combined heuristic divisors to 1A-1B schedule 2025-10-28 22:55:17 +01:00
Marek Lemańczyk 41bff698e8 Add combined heuristic prime tester with shared divisor arrays 2025-10-28 22:40:51 +01:00
mlemanczyk 0850227881 Switched to heuristic GPU for a better performance than the standard one. 2025-10-28 22:15:33 +01:00
Marek Lemańczyk c693a05708 Add reusable group B divisor tables without trailing-three primes 2025-10-28 21:56:18 +01:00
Marek Lemańczyk b83bddf49a Merge pull request #247 from mlemanczyk/codex/fix-unhandled-exception-in-evenperfectbitscanner
Fix GPU Montgomery multiply to avoid unsupported intrinsics
2025-10-28 21:10:09 +01:00
mlemanczyk e5f5cd563e Merge branch 'codex/fix-unhandled-exception-in-evenperfectbitscanner' of https://github.com/mlemanczyk/even-perfect-numbers-scanner into codex/fix-unhandled-exception-in-evenperfectbitscanner 2025-10-28 21:08:14 +01:00
Marek Lemańczyk 3304335072 Rename CPU Montgomery multiply and align GPU helper 2025-10-28 21:06:40 +01:00
mlemanczyk ab045a62c4 Revert "Temporarily enabling heuristic GPU order checks."
This reverts commit 4db06cfe80.
2025-10-28 21:04:09 +01:00
Marek Lemańczyk eaee599f25 Use GPU-safe Montgomery multiply in kernels 2025-10-28 20:58:02 +01:00
mlemanczyk 4db06cfe80 Temporarily enabling heuristic GPU order checks. 2025-10-28 20:49:41 +01:00
mlemanczyk bead40713c Switching .IsPrime check to CPU for much higher performance. 2025-10-28 20:49:11 +01:00
mlemanczyk ca8ab2bcf8 Tester made static to better mimic production code. 2025-10-28 20:48:35 +01:00
Marek Lemańczyk d149700311 Merge pull request #246 from mlemanczyk/codex/compare-cpu-and-gpu-execution-paths
Make GPU sieve results authoritative
2025-10-28 20:14:33 +01:00
Marek Lemańczyk 1792349bad Make GPU sieve results authoritative 2025-10-28 20:14:22 +01:00
mlemanczyk 2f461ee001 Merge branch 'codex/analyze-gc-instances-in-checkdivisors-method' of https://github.com/mlemanczyk/even-perfect-numbers-scanner into codex/analyze-gc-instances-in-checkdivisors-method 2025-10-28 20:12:45 +01:00
mlemanczyk 93aed56f23 Non-heuristic GPU benchmark uncommented, since we improved it. 2025-10-28 20:12:41 +01:00
Marek Lemańczyk ff9774925e Merge pull request #245 from mlemanczyk/codex/move-buildsmallprimes-to-divisorgenerator-g256a9
Enable residue-aware GPU sieve in PrimeTester
2025-10-28 20:00:22 +01:00
Marek Lemańczyk 1f9421430b Reuse GPU divisor tables across PrimeTester batches 2025-10-28 19:53:11 +01:00
Marek Lemańczyk c53966e481 Enable residue-aware GPU sieve in PrimeTester 2025-10-28 18:39:18 +01:00
mlemanczyk 3d43d52f2e Console logging updated 2025-10-28 16:04:42 +01:00
Marek Lemańczyk fa4bb0bcc1 Merge pull request #243 from mlemanczyk/codex/modify-evenperfectbitscanner-to-use-taskthreadpool
Refactor by-divisor scheduling to use task thread pool
2025-10-28 15:47:47 +01:00
Marek Lemańczyk c4c783610b Add start gate to by-divisor task scheduling 2025-10-28 15:32:15 +01:00
Marek Lemańczyk 29f58e3059 Refactor by-divisor scheduling to use task thread pool 2025-10-28 15:11:21 +01:00
Marek Lemańczyk 29af0cb0ce Merge pull request #242 from mlemanczyk/codex/modify-taskthreadpool-to-terminate-threads-after-tasks
Stop tracking double disposal in TaskThreadPool
2025-10-28 14:29:48 +01:00
Marek Lemańczyk 3edf2cd575 Stop tracking double disposal in TaskThreadPool 2025-10-28 14:29:34 +01:00
Marek Lemańczyk a89ed7576e Merge pull request #241 from mlemanczyk/codex/remove-executionscope-structure
Remove GPU context creation lock
2025-10-28 13:39:43 +01:00
Marek Lemańczyk 86990ae0f7 Remove kernel initialization lock from GPU leases 2025-10-28 13:30:47 +01:00
Marek Lemańczyk de7ddd8db3 Remove GPU context creation lock 2025-10-28 13:26:20 +01:00
mlemanczyk 139899fa5e More console logs added. 2025-10-28 02:11:24 +01:00
Marek Lemańczyk c39a99a4fd Merge pull request #240 from mlemanczyk/codex/analyze-mersenne=bydivisor-loops-for-early-breaks-kynlfe
Remove temporary by-divisor review documents
2025-10-28 02:01:52 +01:00
Marek Lemańczyk 6dd5cea79f Remove temporary by-divisor review documents 2025-10-28 02:01:28 +01:00
Marek Lemańczyk bd0727131f Merge pull request #238 from mlemanczyk/codex/optimize-evenperfectbitscanner-loops-kit1so
Hoist prime order thresholds in CPU heuristics
2025-10-27 21:31:23 +01:00
Marek Lemańczyk 968dd89d0d Remove by-divisor path map document 2025-10-27 21:28:41 +01:00
Marek Lemańczyk 94a4ed8328 Hoist prime order thresholds in CPU heuristics 2025-10-27 21:24:43 +01:00
mlemanczyk d381abe0cc We'll now clear the dictionary after renting, to avoid unnecessary clearing if it's not going to be used anymore. 2025-10-27 16:54:49 +01:00
mlemanczyk 50d66a5186 Switched to GPU version in some places. 2025-10-27 15:30:10 +01:00
65 changed files with 5080 additions and 3052 deletions

No files matched your search

@@ -32,7 +32,7 @@ public class GpuUInt128Montgomery64Benchmarks
[Benchmark(Baseline = true)]
public ulong ExtensionMontgomeryMultiply()
{
return Input.Left.MontgomeryMultiply(Input.Right, Input.Modulus, Input.NPrime);
return Input.Left.MontgomeryMultiplyCpu(Input.Right, Input.Modulus, Input.NPrime);
}
/// <summary>
@@ -640,7 +640,7 @@ public class Pow2ModBenchmarks
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
@@ -736,22 +736,22 @@ public class Pow2ModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent);
for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
if (((exponent >> bitIndex) & 1UL) != 0UL)
{
result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
}
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor)
@@ -772,35 +772,35 @@ public class Pow2ModBenchmarks
ulong nibble = remainingExponent & 0xFUL;
ulong base0 = baseVal;
ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime);
ulong base1 = base0.MontgomeryMultiplyCpu(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime);
if ((nibble & 1UL) != 0UL)
{
result = result.MontgomeryMultiply(base0, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime);
}
if ((nibble & 2UL) != 0UL)
{
result = result.MontgomeryMultiply(base1, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime);
}
if ((nibble & 4UL) != 0UL)
{
result = result.MontgomeryMultiply(base2, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime);
}
if ((nibble & 8UL) != 0UL)
{
result = result.MontgomeryMultiply(base3, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base3, modulus, nPrime);
}
baseVal = base3.MontgomeryMultiply(base3, modulus, nPrime);
baseVal = base3.MontgomeryMultiplyCpu(base3, modulus, nPrime);
remainingExponent >>= 4;
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor)
@@ -817,16 +817,16 @@ public class Pow2ModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
const int WindowSize = 5;
Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)];
oddPowers[0] = baseVal;
ulong baseSquared = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
ulong baseSquared = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
for (int i = 1; i < oddPowers.Length; i++)
{
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(baseSquared, modulus, nPrime);
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiplyCpu(baseSquared, modulus, nPrime);
}
int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent);
@@ -835,7 +835,7 @@ public class Pow2ModBenchmarks
{
if (((exponent >> bitIndex) & 1UL) == 0UL)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
bitIndex--;
continue;
}
@@ -852,14 +852,14 @@ public class Pow2ModBenchmarks
int squares = bitIndex - windowStart + 1;
for (int i = 0; i < squares; i++)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
}
result = result.MontgomeryMultiply(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
result = result.MontgomeryMultiplyCpu(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
bitIndex = windowStart - 1;
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
private static ulong Pow2ModBinary(ulong exponent, ulong modulus)
@@ -262,7 +262,7 @@ public class Pow2MontgomeryModBenchmarks
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
@@ -300,22 +300,22 @@ public class Pow2MontgomeryModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent);
for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
if (((exponent >> bitIndex) & 1UL) != 0UL)
{
result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
}
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor)
@@ -336,35 +336,35 @@ public class Pow2MontgomeryModBenchmarks
ulong nibble = remainingExponent & 0xFUL;
ulong base0 = baseVal;
ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime);
ulong base1 = base0.MontgomeryMultiplyCpu(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime);
if ((nibble & 1UL) != 0UL)
{
result = result.MontgomeryMultiply(base0, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime);
}
if ((nibble & 2UL) != 0UL)
{
result = result.MontgomeryMultiply(base1, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime);
}
if ((nibble & 4UL) != 0UL)
{
result = result.MontgomeryMultiply(base2, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime);
}
if ((nibble & 8UL) != 0UL)
{
result = result.MontgomeryMultiply(base3, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(base3, modulus, nPrime);
}
baseVal = base3.MontgomeryMultiply(base3, modulus, nPrime);
baseVal = base3.MontgomeryMultiplyCpu(base3, modulus, nPrime);
remainingExponent >>= 4;
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
// TODO: Update it with the latest implementations
@@ -382,16 +382,16 @@ public class Pow2MontgomeryModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
const int WindowSize = 5;
Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)];
oddPowers[0] = baseVal;
ulong baseSquared = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
ulong baseSquared = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
for (int i = 1; i < oddPowers.Length; i++)
{
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(baseSquared, modulus, nPrime);
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiplyCpu(baseSquared, modulus, nPrime);
}
int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent);
@@ -400,7 +400,7 @@ public class Pow2MontgomeryModBenchmarks
{
if (((exponent >> bitIndex) & 1UL) == 0UL)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
bitIndex--;
continue;
}
@@ -417,14 +417,14 @@ public class Pow2MontgomeryModBenchmarks
int squares = bitIndex - windowStart + 1;
for (int i = 0; i < squares; i++)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
}
result = result.MontgomeryMultiply(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
result = result.MontgomeryMultiplyCpu(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
bitIndex = windowStart - 1;
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
// TODO: Update it with the latest implementations, if it exists in production code. Otherwise just remove TODO
@@ -12,8 +12,8 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser]
public class PrimeTesterBenchmarks
{
private HeuristicPrimeTester _cpuTester = null!;
private HeuristicPrimeTester _gpuTester = null!;
private static readonly HeuristicCombinedPrimeTester _heuristicCombinedTester = new();
private static readonly HeuristicPrimeTester _heuristicTester = new();
private ulong[] _candidates = [];
[ParamsSource(nameof(GetBenchmarkCases))]
@@ -22,8 +22,6 @@ public class PrimeTesterBenchmarks
[GlobalSetup]
public void GlobalSetup()
{
_cpuTester = new HeuristicPrimeTester();
_gpuTester = new HeuristicPrimeTester();
_candidates = BenchmarkCase.Candidates;
if (_candidates.Length == 0)
@@ -31,23 +29,21 @@ public class PrimeTesterBenchmarks
return;
}
CancellationToken cancellationToken = CancellationToken.None;
_cpuTester.IsPrimeCpu(_candidates[0], cancellationToken);
_gpuTester.IsPrimeGpu(_candidates[0]);
_ = _heuristicTester.IsPrimeCpu(_candidates[0]);
_ = _heuristicTester.IsPrimeGpu(_candidates[0]);
_ = Prime.Numbers.IsPrime(_candidates[0]);
}
[Benchmark(Baseline = true)]
public int HeuristicCpu()
{
HeuristicPrimeTester tester = _cpuTester;
HeuristicPrimeTester tester = _heuristicTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (tester.IsPrimeCpu(values[i], cancellationToken))
if (tester.IsPrimeCpu(values[i]))
{
primeCount++;
}
@@ -59,9 +55,8 @@ public class PrimeTesterBenchmarks
[Benchmark]
public int HeuristicGpu()
{
HeuristicPrimeTester tester = _gpuTester;
HeuristicPrimeTester tester = _heuristicTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
@@ -75,18 +70,52 @@ public class PrimeTesterBenchmarks
return primeCount;
}
// This benchmark is very slow. No point in comparing it.
// [Benchmark]
public int NonHeuristicGpu()
[Benchmark]
public int HeuristicCombinedGpu()
{
// PrimeTester tester = _gpuTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (PrimeTester.Exclusive.IsPrimeGpu(values[i], cancellationToken))
if (HeuristicCombinedPrimeTester.IsPrimeGpu(values[i]))
{
primeCount++;
}
}
return primeCount;
}
[Benchmark]
public int HeuristicCombinedCpu()
{
ulong[] values = _candidates;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (HeuristicCombinedPrimeTester.IsPrimeCpu(values[i]))
{
primeCount++;
}
}
return primeCount;
}
private static readonly PrimeTester _nonHeuristicTester = new();
[Benchmark]
public int NonHeuristicGpu()
{
PrimeTester tester = _nonHeuristicTester;
ulong[] values = _candidates;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (tester.IsPrimeGpu(values[i]))
{
primeCount++;
}
@@ -98,14 +127,12 @@ public class PrimeTesterBenchmarks
[Benchmark]
public int NonHeuristicCpu()
{
// PrimeTester tester = _gpuTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (PrimeTester.IsPrimeCpu(values[i], cancellationToken))
if (PrimeTester.IsPrime(values[i]))
{
primeCount++;
}
+1 -1
View File
@@ -55,7 +55,7 @@ public class ProgramTests
[InlineData(81UL, false)]
public void IsPrime_identifies_primes_correctly(ulong n, bool expected)
{
new PrimeTester().IsPrime(n, CancellationToken.None).Should().Be(expected);
PrimeTester.IsPrime(n).Should().Be(expected);
}
[Fact]
+45 -34
View File
@@ -202,14 +202,15 @@ internal readonly struct CliArguments
if (argument.StartsWith("--prime=", StringComparison.OrdinalIgnoreCase))
{
startPrime = Utf8CliParser.ParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1));
ReadOnlySpan<char> primeValue = argument.AsSpan("--prime=".Length);
startPrime = Utf8CliParser.ParseUInt64(primeValue);
startPrimeProvided = true;
continue;
}
if (argument.StartsWith("--max-prime=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
ReadOnlySpan<char> value = argument.AsSpan("--max-prime=".Length);
if (UInt128.TryParse(value, NumberStyles.None, CultureInfo.InvariantCulture, out UInt128 parsedMaxPrime))
{
maxPrimeLimit = parsedMaxPrime;
@@ -238,13 +239,13 @@ internal readonly struct CliArguments
if (argument.StartsWith("--threads=", StringComparison.OrdinalIgnoreCase))
{
threadCount = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
threadCount = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--threads=".Length)));
continue;
}
if (argument.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
ReadOnlySpan<char> value = argument.AsSpan("--mersenne=".Length);
if (value.Equals("pow2mod", StringComparison.OrdinalIgnoreCase))
{
kernelType = GpuKernelType.Pow2Mod;
@@ -281,7 +282,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedLimit))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--divisor-cycles-limit=".Length), out ulong parsedLimit))
{
divisorCyclesSearchLimit = parsedLimit;
}
@@ -291,7 +292,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedResidueMax))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--residue-max-k=".Length), out ulong parsedResidueMax))
{
residueKMax = parsedResidueMax;
}
@@ -301,7 +302,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase))
{
useMersenneOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useMersenneOnGpu = !argument.AsSpan("--mersenne-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
@@ -319,7 +320,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedLimit))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--order-warmup-limit=".Length), out ulong parsedLimit))
{
orderWarmupLimitOverride = parsedLimit;
}
@@ -335,14 +336,14 @@ internal readonly struct CliArguments
if (argument.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
ReadOnlySpan<char> value = argument.AsSpan("--ntt=".Length);
nttBackend = value.Equals("staged", StringComparison.OrdinalIgnoreCase) ? NttBackend.Staged : NttBackend.Reference;
continue;
}
if (argument.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
ReadOnlySpan<char> value = argument.AsSpan("--mod-reduction=".Length);
if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase))
{
modReductionMode = ModReductionMode.Mont64;
@@ -365,25 +366,25 @@ internal readonly struct CliArguments
if (argument.StartsWith("--primes-device=", StringComparison.OrdinalIgnoreCase))
{
forcePrimeKernelsOnCpu = argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
forcePrimeKernelsOnCpu = argument.AsSpan("--primes-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--order-device=", StringComparison.OrdinalIgnoreCase))
{
useOrderOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useOrderOnGpu = !argument.AsSpan("--order-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--rle-blacklist=", StringComparison.OrdinalIgnoreCase))
{
rleBlacklistPath = argument[(argument.IndexOf('=') + 1)..];
rleBlacklistPath = argument["--rle-blacklist=".Length..];
continue;
}
if (argument.StartsWith("--rle-hard-max=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedMax))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--rle-hard-max=".Length), out ulong parsedMax))
{
rleHardMaxP = parsedMax;
}
@@ -393,14 +394,14 @@ internal readonly struct CliArguments
if (argument.StartsWith("--rle-only-last7=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
ReadOnlySpan<char> value = argument.AsSpan("--rle-only-last7=".Length);
rleOnlyLast7 = !value.Equals("false", StringComparison.OrdinalIgnoreCase) && !value.Equals("0", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--zero-hard=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseDouble(argument.AsSpan(argument.IndexOf('=') + 1), out double zeroFraction))
if (Utf8CliParser.TryParseDouble(argument.AsSpan("--zero-hard=".Length), out double zeroFraction))
{
zeroFractionHard = zeroFraction;
}
@@ -410,7 +411,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--zero-conj=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
ReadOnlySpan<char> value = argument.AsSpan("--zero-conj=".Length);
int colon = value.IndexOf(':');
if (colon > 0)
{
@@ -427,87 +428,97 @@ internal readonly struct CliArguments
if (argument.StartsWith("--lucas=", StringComparison.OrdinalIgnoreCase))
{
useMersenneOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useMersenneOnGpu = !argument.AsSpan("--lucas=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase) ||
argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase) ||
argument.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase))
if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase))
{
forcePrimeKernelsOnCpu = argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
forcePrimeKernelsOnCpu = argument.AsSpan("--primes=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase))
{
forcePrimeKernelsOnCpu = argument.AsSpan("--gpu-kernels=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase))
{
forcePrimeKernelsOnCpu = argument.AsSpan("--accelerator=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--results-dir=", StringComparison.OrdinalIgnoreCase))
{
resultsDirectory = argument[(argument.IndexOf('=') + 1)..];
resultsDirectory = argument["--results-dir=".Length..];
continue;
}
if (argument.StartsWith("--results-prefix=", StringComparison.OrdinalIgnoreCase))
{
resultsPrefix = argument[(argument.IndexOf('=') + 1)..];
resultsPrefix = argument["--results-prefix=".Length..];
continue;
}
if (argument.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase))
{
gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-threads=".Length)));
continue;
}
if (argument.StartsWith("--gpu-prime-batch=", StringComparison.OrdinalIgnoreCase))
{
gpuPrimeBatch = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
gpuPrimeBatch = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-batch=".Length)));
continue;
}
if (argument.StartsWith("--ll-slice=", StringComparison.OrdinalIgnoreCase))
{
sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--ll-slice=".Length)));
continue;
}
if (argument.StartsWith("--gpu-scan-batch=", StringComparison.OrdinalIgnoreCase))
{
scanBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
scanBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-scan-batch=".Length)));
continue;
}
if (argument.StartsWith("--block-size=", StringComparison.OrdinalIgnoreCase))
{
blockSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
blockSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--block-size=".Length)));
continue;
}
if (argument.StartsWith("--filter-p=", StringComparison.OrdinalIgnoreCase))
{
filterFile = argument[(argument.IndexOf('=') + 1)..];
filterFile = argument["--filter-p=".Length..];
continue;
}
if (argument.StartsWith("--write-batch-size=", StringComparison.OrdinalIgnoreCase))
{
writeBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
writeBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--write-batch-size=".Length)));
continue;
}
if (argument.StartsWith("--divisor-cycles=", StringComparison.OrdinalIgnoreCase))
{
cyclesPath = argument[(argument.IndexOf('=') + 1)..];
cyclesPath = argument["--divisor-cycles=".Length..];
continue;
}
if (argument.StartsWith("--divisor-cycles-device=", StringComparison.OrdinalIgnoreCase))
{
useGpuCycles = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useGpuCycles = !argument.AsSpan("--divisor-cycles-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--divisor-cycles-batch=", StringComparison.OrdinalIgnoreCase))
{
cyclesBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
cyclesBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--divisor-cycles-batch=".Length)));
continue;
}
@@ -14,7 +14,7 @@ internal static class CalculationResultsFile
using FileStream readStream = new(candidateFile, FileMode.Open, FileAccess.Read, FileShare.ReadWrite);
using StreamReader reader = new(readStream);
string? line;
ReadOnlySpan<char> span;
ReadOnlySpan<char> trimmed;
while ((line = reader.ReadLine()) is not null)
{
if (string.IsNullOrWhiteSpace(line))
@@ -22,36 +22,67 @@ internal static class CalculationResultsFile
continue;
}
span = line.AsSpan().Trim();
if (span.IsEmpty || span[0] == '#')
trimmed = line.AsSpan().Trim();
if (trimmed.IsEmpty || trimmed[0] == '#')
{
continue;
}
int length = trimmed.Length;
int index = 0;
while (index < span.Length)
if (!maxPrimeConfigured)
{
if (!char.IsDigit(span[index]))
while (index < length)
{
index++;
continue;
}
char current = trimmed[index];
uint digit = (uint)(current - '0');
if (digit > 9U)
{
index++;
continue;
}
int start = index;
index++;
while (index < span.Length && char.IsDigit(span[index]))
{
index++;
}
int start = index++;
while (index < length && (uint)(trimmed[index] - '0') <= 9U)
{
index++;
}
if (Utf8CliParser.TryParseUInt64(span[start..index], out ulong parsed))
{
if (!maxPrimeConfigured || (UInt128)parsed <= maxPrimeLimit)
if (Utf8CliParser.TryParseUInt64(trimmed[start..index], out ulong parsed))
{
candidates.Add(parsed);
}
else
}
}
else
{
UInt128 limit = maxPrimeLimit;
while (index < length)
{
char current = trimmed[index];
uint digit = (uint)(current - '0');
if (digit > 9U)
{
skippedByLimit++;
index++;
continue;
}
int start = index++;
while (index < length && (uint)(trimmed[index] - '0') <= 9U)
{
index++;
}
if (Utf8CliParser.TryParseUInt64(trimmed[start..index], out ulong parsed))
{
if ((UInt128)parsed <= limit)
{
candidates.Add(parsed);
}
else
{
skippedByLimit++;
}
}
}
}
+14 -19
View File
@@ -28,6 +28,8 @@ internal static class Program
[ThreadStatic]
private static bool _lastCompositeP;
private static bool _runPrimesOnCpu;
private static void Main(string[] args)
{
try
@@ -49,18 +51,21 @@ internal static class Program
NttGpuMath.GpuTransformBackend = _cliArguments.NttBackend;
NttGpuMath.ReductionMode = _cliArguments.ModReductionMode;
GpuContextPool.ForceCpu = _cliArguments.ForcePrimeKernelsOnCpu;
_runPrimesOnCpu = _cliArguments.ForcePrimeKernelsOnCpu;
ulong currentP = _cliArguments.StartPrime;
ulong remainder = currentP % 6UL;
bool startPrimeProvided = _cliArguments.StartPrimeProvided;
int threadCount = Math.Max(1, _cliArguments.ThreadCount);
int gpuPrimeBatch = Math.Max(1, _cliArguments.GpuPrimeBatch);
UnboundedTaskScheduler.ConfigureThreadCount(threadCount);
PrimeTester.GpuBatchSize = gpuPrimeBatch;
GpuContextPool.WarmUpPool(threadCount);
PrimeTester.WarmUpGpuKernels(threadCount);
Console.WriteLine("Starting up threads...");
_ = UnboundedTaskScheduler.Instance;
int blockSize = Math.Max(1, _cliArguments.BlockSize);
int gpuPrimeThreads = Math.Max(1, _cliArguments.GpuPrimeThreads);
int gpuPrimeBatch = Math.Max(1, _cliArguments.GpuPrimeBatch);
GpuKernelType kernelType = _cliArguments.KernelType;
bool useBitTransform = _cliArguments.UseBitTransform;
bool useOrder = _cliArguments.UseOrder;
@@ -104,14 +109,7 @@ internal static class Program
testPrimeCandidateLimit = Math.Max(1, threadCount) * 3;
}
if (useByDivisor)
{
_byDivisorStartPrime = startPrimeProvided ? currentP : 0UL;
}
else
{
_byDivisorStartPrime = 0UL;
}
_byDivisorStartPrime = useByDivisor && startPrimeProvided ? currentP : 0UL;
if (useByDivisor && !testMode && string.IsNullOrEmpty(filterFile))
{
@@ -121,7 +119,7 @@ internal static class Program
// Apply GPU prime sieve runtime configuration
GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads);
PrimeTester.GpuBatchSize = Math.Max(1, gpuPrimeBatch);
PrimeTester.GpuBatchSize = gpuPrimeBatch;
MersenneDivisorCycles mersenneDivisorCycles = new();
if (!File.Exists(cyclesPath) || continueCyclesGeneration)
@@ -177,10 +175,7 @@ internal static class Program
Console.WriteLine("Divisor cycles are ready");
if (useByDivisor)
{
blockSize = 1;
}
blockSize = useByDivisor ? 1 : blockSize;
CandidatesCalculator.Initialize(currentP);
PrimeTesters = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
@@ -266,7 +261,7 @@ internal static class Program
orderWarmupLimitOverride ?? 5_000_000UL,
NttGpuMath.ReductionMode,
mersenneOnGpu ? "gpu" : "cpu",
GpuContextPool.ForceCpu ? "cpu" : "gpu",
_runPrimesOnCpu ? "cpu" : "gpu",
orderOnGpu ? "gpu" : "cpu");
if (!string.IsNullOrEmpty(_cliArguments.ResultsPrefix))
@@ -631,9 +626,9 @@ internal static class Program
// If primes-device=gpu, route p primality through GPU-assisted sieve with deterministic MR validation.
if (!_cliArguments.UseByDivisor)
{
if (!(GpuContextPool.ForceCpu
? PrimeTesters.Value!.IsPrime(p, CancellationToken.None)
: PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None)))
if (!(_runPrimesOnCpu
? PrimeTester.IsPrime(p)
: PrimeTesters.Value!.IsPrimeGpu(p)))
{
_lastCompositeP = true;
return false;
@@ -9,88 +9,88 @@ namespace PerfectNumbers.Core.Tests.Gpu;
[Collection("GpuNtt")]
public class GpuContextPoolTests
{
[Fact]
[Trait("Category", "Fast")]
public void RentPreferred_cpu_returns_cpu_accelerator()
{
GpuContextPool.DisposeAll();
// [Fact]
// [Trait("Category", "Fast")]
// public void RentPreferred_cpu_returns_cpu_accelerator()
// {
// GpuContextPool.DisposeAll();
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
}
finally
{
lease.Dispose();
}
// GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
// }
// finally
// {
// lease.Dispose();
// }
GpuContextPool.DisposeAll();
}
// GpuContextPool.DisposeAll();
// }
[Fact]
[Trait("Category", "Fast")]
public void RentPreferred_reuses_cpu_context()
{
GpuContextPool.DisposeAll();
// [Fact]
// [Trait("Category", "Fast")]
// public void RentPreferred_reuses_cpu_context()
// {
// GpuContextPool.DisposeAll();
Context firstContext;
Accelerator firstAccelerator;
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
firstContext = lease.Context;
firstAccelerator = lease.Accelerator;
}
finally
{
lease.Dispose();
}
// Context firstContext;
// Accelerator firstAccelerator;
// GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// firstContext = lease.Context;
// firstAccelerator = lease.Accelerator;
// }
// finally
// {
// lease.Dispose();
// }
lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
lease.Context.Should().BeSameAs(firstContext);
lease.Accelerator.Should().BeSameAs(firstAccelerator);
}
finally
{
lease.Dispose();
}
// lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// lease.Context.Should().BeSameAs(firstContext);
// lease.Accelerator.Should().BeSameAs(firstAccelerator);
// }
// finally
// {
// lease.Dispose();
// }
GpuContextPool.DisposeAll();
}
// GpuContextPool.DisposeAll();
// }
[Fact]
[Trait("Category", "Fast")]
public void DisposeAll_clears_pooled_contexts()
{
GpuContextPool.DisposeAll();
// [Fact]
// [Trait("Category", "Fast")]
// public void DisposeAll_clears_pooled_contexts()
// {
// GpuContextPool.DisposeAll();
Context firstContext;
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
firstContext = lease.Context;
}
finally
{
lease.Dispose();
}
// Context firstContext;
// GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// firstContext = lease.Context;
// }
// finally
// {
// lease.Dispose();
// }
GpuContextPool.DisposeAll();
// GpuContextPool.DisposeAll();
lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
lease.Context.Should().NotBeSameAs(firstContext);
}
finally
{
lease.Dispose();
}
// lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// lease.Context.Should().NotBeSameAs(firstContext);
// }
// finally
// {
// lease.Dispose();
// }
GpuContextPool.DisposeAll();
}
// GpuContextPool.DisposeAll();
// }
}
@@ -121,7 +121,7 @@ public class GpuKernelCompilationTests
private static readonly Action<Index1D, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> MersenneDivisorCyclesKernel =
DivisorCycleKernels.GpuDivisorCycleKernel;
private static readonly Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>> PrimeTesterSmallPrimeKernel =
private static readonly Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> PrimeTesterSmallPrimeKernel =
PrimeTesterKernels.SmallPrimeSieveKernel;
private static readonly Action<Index1D, ArrayView<ulong>, ArrayView<byte>> PrimeTesterSharesFactorKernel =
@@ -167,23 +167,23 @@ public class GpuKernelCompilationTests
return new object[] { name, loader };
}
[Theory]
[Trait("Category", "Fast")]
[MemberData(nameof(KernelLoaders))]
public void Kernel_compiles_without_internal_compiler_errors(string name, Action<Accelerator> compile)
{
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
var accelerator = lease.Accelerator;
Action action = () => compile(accelerator);
action.Should().NotThrow($"kernel {name} should compile");
}
finally
{
lease.Dispose();
}
}
// [Theory]
// [Trait("Category", "Fast")]
// [MemberData(nameof(KernelLoaders))]
// public void Kernel_compiles_without_internal_compiler_errors(string name, Action<Accelerator> compile)
// {
// GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// var accelerator = lease.Accelerator;
// Action action = () => compile(accelerator);
// action.Should().NotThrow($"kernel {name} should compile");
// }
// finally
// {
// lease.Dispose();
// }
// }
private static void CompileDivisorCycleCacheKernel(Accelerator accelerator)
{
@@ -336,7 +336,7 @@ public class GpuKernelCompilationTests
private static void CompilePrimeTesterSmallPrimeKernel(Accelerator accelerator)
{
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>>(PrimeTesterSmallPrimeKernel);
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterSmallPrimeKernel);
_ = kernel;
}
@@ -18,22 +18,33 @@ public class GpuPrimeLimiterTests
var startedSecond = new TaskCompletionSource<bool>();
var sw = new Stopwatch();
using var first = GpuPrimeWorkLimiter.Acquire();
var first = GpuPrimeWorkLimiter.Acquire();
sw.Start();
var secondTask = Task.Run(() =>
{
using var second = GpuPrimeWorkLimiter.Acquire();
startedSecond.TrySetResult(true);
var second = GpuPrimeWorkLimiter.Acquire();
try
{
startedSecond.TrySetResult(true);
}
finally
{
second.Dispose();
}
});
// Give secondTask a moment to attempt acquire (it should block).
await Task.Delay(50);
startedSecond.Task.IsCompleted.Should().BeFalse();
// Release first after short delay, allowing second to continue.
// (Dispose at scope end, but we want to end now.)
first.Dispose();
try
{
// Give secondTask a moment to attempt acquire (it should block).
await Task.Delay(50);
startedSecond.Task.IsCompleted.Should().BeFalse();
}
finally
{
// Release first after the initial wait so the second task can proceed.
first.Dispose();
}
await secondTask.WaitAsync(TimeSpan.FromSeconds(1));
sw.Stop();
@@ -14,7 +14,7 @@ public class LucasLehmerPrimeTesterTests
[Trait("Category", "Fast")]
public void IsMersennePrime_gpu_returns_expected_results(ulong p, bool expected)
{
new MersenneNumberLucasLehmerGpuTester().IsMersennePrime(p).Should().Be(expected);
new MersenneNumberLucasLehmerGpuTester().IsPrime(p, true).Should().Be(expected);
}
[Fact]
@@ -26,7 +26,7 @@ public class LucasLehmerPrimeTesterTests
const ulong exponent = 128UL;
var tester = new MersenneNumberLucasLehmerGpuTester();
tester.WarmUpNttParameters(exponent);
tester.IsMersennePrime(exponent).Should().BeFalse();
tester.IsPrime(exponent, true).Should().BeFalse();
}
[Fact]
@@ -49,7 +49,7 @@ public class LucasLehmerPrimeTesterTests
{
var tester = new MersenneNumberLucasLehmerGpuTester();
tester.WarmUpNttParameters(exponent);
tester.IsMersennePrime(exponent).Should().BeFalse();
tester.IsPrime(exponent, true).Should().BeFalse();
}
[Theory]
@@ -10,19 +10,17 @@ public class PrimeTesterGpuTests
[Fact]
public void IsPrimeGpu_accepts_known_primes()
{
GpuContextPool.ForceCpu = false;
try
{
var tester = new PrimeTester();
ulong[] primes = [31UL, 61UL, 89UL, 107UL, 127UL, 521UL];
foreach (ulong prime in primes)
{
tester.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue();
tester.IsPrimeGpu(prime).Should().BeTrue();
}
}
finally
{
GpuContextPool.ForceCpu = false;
GpuContextPool.DisposeAll();
}
}
@@ -30,19 +28,17 @@ public class PrimeTesterGpuTests
[Fact]
public void IsPrimeGpu_rejects_composites()
{
GpuContextPool.ForceCpu = false;
try
{
var tester = new PrimeTester();
ulong[] composites = [33UL, 39UL, 51UL, 77UL, 91UL];
foreach (ulong composite in composites)
{
tester.IsPrimeGpu(composite, CancellationToken.None).Should().BeFalse();
tester.IsPrimeGpu(composite).Should().BeFalse();
}
}
finally
{
GpuContextPool.ForceCpu = false;
GpuContextPool.DisposeAll();
}
}
+18 -18
View File
@@ -7,26 +7,26 @@ namespace PerfectNumbers.Core.Tests;
[Collection("GpuNtt")]
public class GpuKernelPoolTests
{
[Fact]
[Trait("Category", "Fast")]
public void EnsureSmallCyclesOnDevice_returns_same_view_on_subsequent_calls()
{
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
var accelerator = lease.Accelerator;
// [Fact]
// [Trait("Category", "Fast")]
// public void EnsureSmallCyclesOnDevice_returns_same_view_on_subsequent_calls()
// {
// GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
// try
// {
// var accelerator = lease.Accelerator;
var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
var second = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
// var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
// var second = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
second.Equals(first).Should().BeTrue();
second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
}
finally
{
lease.Dispose();
}
}
// second.Equals(first).Should().BeTrue();
// second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
// }
// finally
// {
// lease.Dispose();
// }
// }
[Fact]
[Trait("Category", "Fast")]
@@ -17,19 +17,31 @@ public class GpuWorkLimiterTests
var startedSecond = new TaskCompletionSource<bool>();
var sw = new Stopwatch();
using var first = GpuWorkLimiter.Acquire();
var first = GpuWorkLimiter.Acquire();
sw.Start();
var secondTask = Task.Run(() =>
{
using var second = GpuWorkLimiter.Acquire();
startedSecond.TrySetResult(true);
var second = GpuWorkLimiter.Acquire();
try
{
startedSecond.TrySetResult(true);
}
finally
{
second.Dispose();
}
});
await Task.Delay(50);
startedSecond.Task.IsCompleted.Should().BeFalse();
first.Dispose();
try
{
await Task.Delay(50);
startedSecond.Task.IsCompleted.Should().BeFalse();
}
finally
{
first.Dispose();
}
await secondTask.WaitAsync(TimeSpan.FromSeconds(1));
sw.Stop();
@@ -0,0 +1,239 @@
using System;
using System.Numerics;
using System.Collections.Generic;
using PerfectNumbers.Core;
using FluentAssertions;
using Xunit;
namespace PerfectNumbers.Core.Tests;
public class HeuristicCombinedPrimeTesterTests
{
[Theory]
[Trait("Category", "Fast")]
[InlineData((byte)1)]
[InlineData((byte)7)]
[InlineData((byte)9)]
public void Enumerator_uses_one_a_one_b_schedule_by_default(byte lastDigit)
{
Span<HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState> groupBState =
stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4];
var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500UL * 500UL, lastDigit, groupBState);
var results = new List<HeuristicCombinedPrimeTester.HeuristicDivisorCandidate>(32);
while (results.Count < 20 && enumerator.TryGetNext(out var candidate))
{
results.Add(candidate);
}
const int segmentsToCheck = 6;
results.Count.Should().BeGreaterThanOrEqualTo(4 + segmentsToCheck * 2);
results[0].Value.Should().Be(3UL);
results[0].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant);
results[1].Value.Should().Be(7UL);
results[1].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant);
results[2].Value.Should().Be(11UL);
results[2].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant);
results[3].Value.Should().Be(13UL);
results[3].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant);
int index = 4;
for (int segment = 0; segment < segmentsToCheck; segment++)
{
HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupACandidate = results[index++];
groupACandidate.Value.Should().BeGreaterThan(13UL);
(groupACandidate.Value % 10UL).Should().Be(3UL);
groupACandidate.Group.Should().NotBe(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB);
HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++];
groupBCandidate.Value.Should().BeGreaterThan(13UL);
(groupBCandidate.Value % 10UL).Should().NotBe(3UL);
groupBCandidate.Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB);
}
}
[Theory]
[Trait("Category", "Fast")]
[InlineData((byte)1)]
[InlineData((byte)7)]
[InlineData((byte)9)]
public void One_a_one_b_pattern_interleaves_group_a_and_group_b_values(byte lastDigit)
{
ReadOnlySpan<ulong> combined = HeuristicCombinedPrimeTester.GetCombinedDivisors(
lastDigit,
HeuristicCombinedPrimeTester.CombinedDivisorPattern.OneAOneB);
combined.Length.Should().BeGreaterThan(18);
combined[0].Should().Be(3UL);
combined[1].Should().Be(7UL);
combined[2].Should().Be(11UL);
combined[3].Should().Be(13UL);
int index = 4;
const int segmentsToCheck = 6;
for (int segment = 0; segment < segmentsToCheck; segment++)
{
ulong groupAValue = combined[index++];
groupAValue.Should().BeGreaterThan(13UL);
(groupAValue % 10UL).Should().Be(3UL);
ulong groupBValue = combined[index++];
groupBValue.Should().BeGreaterThan(13UL);
(groupBValue % 10UL).Should().NotBe(3UL);
}
}
[Theory]
[Trait("Category", "Fast")]
[InlineData((byte)1)]
[InlineData((byte)7)]
[InlineData((byte)9)]
public void Two_a_one_b_pattern_interleaves_two_group_a_values_before_group_b(byte lastDigit)
{
ReadOnlySpan<ulong> combined = HeuristicCombinedPrimeTester.GetCombinedDivisors(
lastDigit,
HeuristicCombinedPrimeTester.CombinedDivisorPattern.TwoAOneB);
combined.Length.Should().BeGreaterThan(20);
combined[0].Should().Be(3UL);
combined[1].Should().Be(7UL);
combined[2].Should().Be(11UL);
combined[3].Should().Be(13UL);
int index = 4;
const int segmentsToCheck = 5;
for (int segment = 0; segment < segmentsToCheck; segment++)
{
for (int i = 0; i < 2; i++)
{
ulong groupAValue = combined[index++];
groupAValue.Should().BeGreaterThan(13UL);
(groupAValue % 10UL).Should().Be(3UL);
}
ulong groupBValue = combined[index++];
groupBValue.Should().BeGreaterThan(13UL);
(groupBValue % 10UL).Should().NotBe(3UL);
}
}
[Theory]
[Trait("Category", "Fast")]
[InlineData((byte)1)]
[InlineData((byte)7)]
[InlineData((byte)9)]
public void Three_a_two_b_pattern_interleaves_three_group_a_values_before_two_group_b_values(byte lastDigit)
{
ReadOnlySpan<ulong> combined = HeuristicCombinedPrimeTester.GetCombinedDivisors(
lastDigit,
HeuristicCombinedPrimeTester.CombinedDivisorPattern.ThreeATwoB);
combined.Length.Should().BeGreaterThan(28);
combined[0].Should().Be(3UL);
combined[1].Should().Be(7UL);
combined[2].Should().Be(11UL);
combined[3].Should().Be(13UL);
int index = 4;
const int segmentsToCheck = 4;
for (int segment = 0; segment < segmentsToCheck; segment++)
{
for (int i = 0; i < 3; i++)
{
ulong groupAValue = combined[index++];
groupAValue.Should().BeGreaterThan(13UL);
(groupAValue % 10UL).Should().Be(3UL);
}
for (int i = 0; i < 2; i++)
{
ulong groupBValue = combined[index++];
groupBValue.Should().BeGreaterThan(13UL);
(groupBValue % 10UL).Should().NotBe(3UL);
}
}
}
[Theory]
[Trait("Category", "Fast")]
[InlineData((byte)1)]
[InlineData((byte)7)]
[InlineData((byte)9)]
public void Three_a_one_b_pattern_interleaves_three_group_a_values_before_group_b(byte lastDigit)
{
ReadOnlySpan<ulong> combined = HeuristicCombinedPrimeTester.GetCombinedDivisors(
lastDigit,
HeuristicCombinedPrimeTester.CombinedDivisorPattern.ThreeAOneB);
combined.Length.Should().BeGreaterThan(24);
combined[0].Should().Be(3UL);
combined[1].Should().Be(7UL);
combined[2].Should().Be(11UL);
combined[3].Should().Be(13UL);
int index = 4;
const int segmentsToCheck = 4;
for (int segment = 0; segment < segmentsToCheck; segment++)
{
for (int i = 0; i < 3; i++)
{
ulong groupAValue = combined[index++];
groupAValue.Should().BeGreaterThan(13UL);
(groupAValue % 10UL).Should().Be(3UL);
}
ulong groupBValue = combined[index++];
groupBValue.Should().BeGreaterThan(13UL);
(groupBValue % 10UL).Should().NotBe(3UL);
}
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(1UL)]
[InlineData(2UL)]
[InlineData(3UL)]
[InlineData(4UL)]
[InlineData(5UL)]
[InlineData(8UL)]
[InlineData(9UL)]
[InlineData(15UL)]
[InlineData(16UL)]
[InlineData(17UL)]
[InlineData(24UL)]
[InlineData(25UL)]
[InlineData(26UL)]
[InlineData(255UL)]
[InlineData(256UL)]
[InlineData(257UL)]
[InlineData(1_000UL)]
[InlineData(1_048_575UL)]
[InlineData(1_048_576UL)]
[InlineData(1_048_577UL)]
[InlineData(4_000_000UL)]
[InlineData(4_000_001UL)]
[InlineData(4_000_004UL)]
[InlineData(65_535_999UL)]
[InlineData(65_536_000UL)]
[InlineData(65_536_001UL)]
[InlineData(4_294_967_295UL)]
[InlineData(4_294_967_296UL)]
[InlineData(4_294_967_297UL)]
[InlineData(18_446_744_073_709_551_615UL)]
[InlineData(18_446_744_073_709_551_614UL)]
[InlineData(18_446_744_073_709_551_613UL)]
public void ComputeHeuristicDivisorSquareLimit_returns_input_value(ulong value)
{
ulong limit = HeuristicCombinedPrimeTester.ComputeHeuristicDivisorSquareLimit(value);
limit.Should().Be(value);
}
}
@@ -30,18 +30,25 @@ public class MersenneNumberDivisorCpuTesterTests
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(13UL);
using var session = tester.CreateDivisorSession();
ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
byte[] hits = new byte[primes.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
byte[] hits = new byte[primes.Length];
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
Array.Fill(hits, (byte)0);
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
Array.Fill(hits, (byte)0);
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -51,15 +58,22 @@ public class MersenneNumberDivisorCpuTesterTests
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(43UL);
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -69,15 +83,22 @@ public class MersenneNumberDivisorCpuTesterTests
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(61UL);
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
finally
{
session.Dispose();
}
}
}
@@ -149,13 +149,20 @@ public class MersenneNumberDivisorGpuTesterTests
var tester = new MersenneNumberDivisorByDivisorGpuTester();
tester.ConfigureFromMaxPrime(43UL);
using var session = tester.CreateDivisorSession();
ulong[] primes = { 31UL, 37UL, 41UL };
byte[] hits = new byte[primes.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL };
byte[] hits = new byte[primes.Length];
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes, hits);
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0 });
hits.Should().ContainInOrder(new byte[] { 0, 1, 0 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -165,18 +172,25 @@ public class MersenneNumberDivisorGpuTesterTests
var tester = new MersenneNumberDivisorByDivisorGpuTester();
tester.ConfigureFromMaxPrime(43UL);
using var session = tester.CreateDivisorSession();
ulong[] primes = { 31UL, 37UL, 41UL, 43UL };
byte[] hits = new byte[primes.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL, 43UL };
byte[] hits = new byte[primes.Length];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 });
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 });
Array.Fill(hits, (byte)0);
ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL));
session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
Array.Fill(hits, (byte)0);
ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL));
session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -190,14 +204,21 @@ public class MersenneNumberDivisorGpuTesterTests
tester.ConfigureFromMaxPrime(47UL);
using var session = tester.CreateDivisorSession();
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -210,16 +231,23 @@ public class MersenneNumberDivisorGpuTesterTests
};
tester.ConfigureFromMaxPrime(43UL);
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -232,16 +260,23 @@ public class MersenneNumberDivisorGpuTesterTests
};
tester.ConfigureFromMaxPrime(61UL);
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
finally
{
session.Dispose();
}
}
[Fact]
@@ -297,14 +332,21 @@ public class MersenneNumberDivisorGpuTesterTests
tester.ConfigureFromMaxPrime(47UL);
using var session = tester.CreateDivisorSession();
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
}
finally
{
session.Dispose();
}
}
}
@@ -17,12 +17,29 @@ public class PrimeTesterIsPrimeTests
[InlineData(97UL, true)]
[InlineData(133UL, false)]
[InlineData(137UL, true)]
[InlineData(143UL, false)]
[InlineData(199UL, true)]
[InlineData(209UL, false)]
[InlineData(403UL, false)]
[InlineData(341UL, false)]
public void IsPrime_returns_expected_results(ulong n, bool expected)
{
var tester = new PrimeTester();
tester.IsPrime(n, CancellationToken.None).Should().Be(expected);
PrimeTester.IsPrime(n).Should().Be(expected);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(139UL, true)]
[InlineData(189UL, false)]
[InlineData(193UL, true)]
[InlineData(279UL, false)]
public void IsPrimeGpu_uses_residue_specific_tables(ulong n, bool expected)
{
var tester = new PrimeTester();
tester.IsPrimeGpu(n).Should().Be(expected);
}
[Fact]
@@ -30,20 +47,14 @@ public class PrimeTesterIsPrimeTests
public void IsPrimeGpu_falls_back_to_cpu_when_forced()
{
var tester = new PrimeTester();
GpuContextPool.ForceCpu = true;
tester.IsPrimeGpu(137UL, CancellationToken.None).Should().BeTrue();
tester.IsPrimeGpu(341UL, CancellationToken.None).Should().BeFalse();
GpuContextPool.ForceCpu = false;
PrimeTester.IsPrime(137UL).Should().BeTrue();
PrimeTester.IsPrime(341UL).Should().BeFalse();
}
[Fact]
[Trait("Category", "Fast")]
public void IsPrimeGpu_handles_parallel_calls_for_known_primes()
{
GpuContextPool.ForceCpu = false;
GpuPrimeWorkLimiter.SetLimit(64);
PrimeTester.GpuBatchSize = 256_000;
@@ -53,7 +64,7 @@ public class PrimeTesterIsPrimeTests
var testers = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
Parallel.For(0, primes.Length, i =>
{
results[i] = testers.Value!.IsPrimeGpu(primes[i], CancellationToken.None);
results[i] = testers.Value!.IsPrimeGpu(primes[i]);
});
results.Should().AllBeEquivalentTo(true);
@@ -62,7 +73,7 @@ public class PrimeTesterIsPrimeTests
{
foreach (ulong prime in primes)
{
testers.Value!.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue();
testers.Value!.IsPrimeGpu(prime).Should().BeTrue();
}
});
}
@@ -32,14 +32,39 @@ public class PrimesGeneratorTests
var lastOnePow2 = PrimesGenerator.SmallPrimesPow2LastOne;
var lastSeven = PrimesGenerator.SmallPrimesLastSeven;
var lastSevenPow2 = PrimesGenerator.SmallPrimesPow2LastSeven;
var lastThree = DivisorGenerator.SmallPrimesLastThree;
var lastThreePow2 = DivisorGenerator.SmallPrimesPow2LastThree;
var lastNine = DivisorGenerator.SmallPrimesLastNine;
var lastNinePow2 = DivisorGenerator.SmallPrimesPow2LastNine;
var lastOneWithoutThree = DivisorGenerator.SmallPrimesLastOneWithoutLastThree;
var lastOneWithoutThreePow2 = DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree;
var lastSevenWithoutThree = DivisorGenerator.SmallPrimesLastSevenWithoutLastThree;
var lastSevenWithoutThreePow2 = DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree;
var lastNineWithoutThree = DivisorGenerator.SmallPrimesLastNineWithoutLastThree;
var lastNineWithoutThreePow2 = DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree;
lastOne.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit);
lastSeven.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit);
lastThree.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit);
lastNine.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit);
lastOneWithoutThree.Length.Should().BeGreaterThan(0);
lastSevenWithoutThree.Length.Should().BeGreaterThan(0);
lastNineWithoutThree.Length.Should().BeGreaterThan(0);
lastOneWithoutThree.Length.Should().BeLessThan(lastOne.Length);
lastSevenWithoutThree.Length.Should().BeLessThan(lastSeven.Length);
lastNineWithoutThree.Length.Should().BeLessThan(lastNine.Length);
lastOnePow2.Length.Should().Be(lastOne.Length);
lastSevenPow2.Length.Should().Be(lastSeven.Length);
lastThreePow2.Length.Should().Be(lastThree.Length);
lastNinePow2.Length.Should().Be(lastNine.Length);
lastOneWithoutThreePow2.Length.Should().Be(lastOneWithoutThree.Length);
lastSevenWithoutThreePow2.Length.Should().Be(lastSevenWithoutThree.Length);
lastNineWithoutThreePow2.Length.Should().Be(lastNineWithoutThree.Length);
lastOne[..5].Should().Equal([3U, 7U, 11U, 13U, 19U]);
lastSeven[..5].Should().Equal([3U, 7U, 11U, 13U, 17U]);
lastThree[..6].Should().Equal([3U, 7U, 11U, 13U, 17U, 19U]);
lastNine[..6].Should().Equal([3U, 7U, 11U, 13U, 17U, 19U]);
for (int i = 0; i < lastOne.Length; i++)
{
@@ -58,6 +83,64 @@ public class PrimesGeneratorTests
.Should().BeTrue();
lastSevenPow2[i].Should().Be(prime * (ulong)prime);
}
for (int i = 0; i < lastThree.Length; i++)
{
uint prime = lastThree[i];
uint mod10 = prime % 10U;
(mod10 == 3U || mod10 == 7U || prime == 11U || prime == 19U)
.Should().BeTrue();
lastThreePow2[i].Should().Be(prime * (ulong)prime);
}
lastOneWithoutThree.Should().NotBeEmpty();
lastSevenWithoutThree.Should().NotBeEmpty();
lastNineWithoutThree.Should().NotBeEmpty();
for (int i = 0; i < lastOneWithoutThree.Length; i++)
{
uint prime = lastOneWithoutThree[i];
(prime % 10U != 3U).Should().BeTrue();
if (i > 0)
{
prime.Should().BeGreaterThan(lastOneWithoutThree[i - 1]);
}
lastOneWithoutThreePow2[i].Should().Be(prime * (ulong)prime);
}
for (int i = 0; i < lastSevenWithoutThree.Length; i++)
{
uint prime = lastSevenWithoutThree[i];
(prime % 10U != 3U).Should().BeTrue();
if (i > 0)
{
prime.Should().BeGreaterThan(lastSevenWithoutThree[i - 1]);
}
lastSevenWithoutThreePow2[i].Should().Be(prime * (ulong)prime);
}
for (int i = 0; i < lastNineWithoutThree.Length; i++)
{
uint prime = lastNineWithoutThree[i];
(prime % 10U != 3U).Should().BeTrue();
if (i > 0)
{
prime.Should().BeGreaterThan(lastNineWithoutThree[i - 1]);
}
lastNineWithoutThreePow2[i].Should().Be(prime * (ulong)prime);
}
for (int i = 0; i < lastNine.Length; i++)
{
uint prime = lastNine[i];
uint mod10 = prime % 10U;
(mod10 == 3U || mod10 == 7U || mod10 == 9U || prime == 11U)
.Should().BeTrue();
lastNinePow2[i].Should().Be(prime * (ulong)prime);
}
}
}
@@ -2,11 +2,8 @@ namespace PerfectNumbers.Core.Cpu;
internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{
private bool _disposed;
public void Reset()
{
_disposed = false;
}
public void CheckDivisor(
@@ -67,12 +64,6 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi
public void Dispose()
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(MersenneCpuDivisorScanSession), "Divisor scan session disposed twice.");
}
_disposed = true;
ThreadStaticPools.ReturnMersenneCpuDivisorSession(this);
}
}
@@ -7,10 +7,7 @@ namespace PerfectNumbers.Core.Cpu;
public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester
{
private readonly object _sync = new();
private ulong _divisorLimit;
private ulong _lastStatusDivisor;
private int _batchSize = 1_024;
[ThreadStatic]
@@ -32,7 +29,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
public void ConfigureFromMaxPrime(ulong maxPrime)
{
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_lastStatusDivisor = 0UL;
}
public ulong DivisorLimit
@@ -60,22 +56,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
// return true;
// }
ulong processedCount;
bool processedAll;
bool composite = CheckDivisors(
prime,
allowedMax,
out processedAll,
out processedCount);
if (processedCount > 0UL)
{
lock (_sync)
{
UpdateStatusUnsafe(processedCount);
}
}
out processedAll);
if (composite)
{
@@ -111,10 +97,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private bool CheckDivisors(
ulong prime,
ulong allowedMax,
out bool processedAll,
out ulong processedCount)
out bool processedAll)
{
processedCount = 0UL;
processedAll = false;
// The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs.
@@ -151,44 +135,78 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
ulong stepHigh = step.High;
ulong stepLow = step.Low;
byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
byte step8 = (byte)(stepLow % 8UL);
byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL);
byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL);
byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL);
byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL);
ulong divisorHigh = divisor.High;
ulong limitHigh = limit.High;
ulong divisorLow = divisor.Low;
byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL);
byte remainder8 = (byte)(divisorLow % 8UL);
byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL);
byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL);
byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL);
// 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;
bool lastIsSeven = lastDigit == LastDigit.Seven;
const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);
byte step10;
byte step8;
byte step3;
byte step7;
byte step11;
byte remainder10;
byte remainder8;
byte remainder3;
byte remainder7;
byte remainder11;
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);
remainder10 = (byte)(divisorLow % 10UL);
remainder8 = (byte)(divisorLow & 7UL);
remainder3 = (byte)(divisorLow % 3UL);
remainder7 = (byte)(divisorLow % 7UL);
remainder11 = (byte)(divisorLow % 11UL);
return CheckDivisors64(
prime,
stepLow,
limit.Low,
divisorLow,
decimalMask,
step10,
step8,
step3,
step7,
step11,
remainder10,
remainder8,
remainder3,
remainder7,
remainder11,
out processedAll);
}
ulong divisorHigh = divisor.High;
step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
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);
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);
while (divisor.CompareTo(limit) <= 0)
{
ulong candidate = divisor.Low;
processedCount++;
ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
bool admissible = ((decimalMask >> remainder10) & 1) != 0;
if (admissible
&& (remainder8 == 1 || remainder8 == 7)
&& remainder3 != 0
&& remainder5 != 0
&& remainder7 != 0
&& remainder11 != 0)
bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0;
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
{
ulong candidate = divisor.Low;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
ulong divisorCycle;
// Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs.
@@ -233,16 +251,137 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
divisor.Add(step);
remainder10 = AddMod(remainder10, step10, (byte)10);
remainder8 = AddMod(remainder8, step8, (byte)8);
remainder5 = AddMod(remainder5, step5, (byte)5);
remainder3 = AddMod(remainder3, step3, (byte)3);
remainder7 = AddMod(remainder7, step7, (byte)7);
remainder11 = AddMod(remainder11, step11, (byte)11);
remainder10 = AddMod10(remainder10, step10);
remainder8 = AddMod8(remainder8, step8);
remainder3 = AddMod3(remainder3, step3);
remainder7 = AddMod7(remainder7, step7);
remainder11 = AddMod11(remainder11, step11);
}
processedAll = true;
return false;
}
private static bool CheckDivisors64(
ulong prime,
ulong step,
ulong limit,
ulong divisor,
ushort decimalMask,
byte step10,
byte step8,
byte step3,
byte step7,
byte step11,
byte remainder10,
byte remainder8,
byte remainder3,
byte remainder7,
byte remainder11,
out bool processedAll)
{
processedAll = false;
bool canAdvance = step <= limit;
if (!canAdvance)
{
if (divisor <= limit)
{
if (remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
{
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
ulong divisorCycle;
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
divisor,
prime,
divisorData,
out ulong computedCycle,
out bool primeOrderFailed) || computedCycle == 0UL)
{
divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
divisor,
divisorData,
skipPrimeOrderHeuristic: primeOrderFailed);
}
else
{
divisorCycle = computedCycle;
}
if (divisorCycle == prime)
{
processedAll = true;
return true;
}
if (divisorCycle == 0UL)
{
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
}
}
}
processedAll = true;
return false;
}
processedAll = divisor.CompareTo(limit) > 0;
return false;
if (divisor > limit)
{
processedAll = true;
return false;
}
ulong remainingIterations = ((limit - divisor) / step) + 1UL;
while (true)
{
bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0;
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
{
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
ulong divisorCycle;
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
divisor,
prime,
divisorData,
out ulong computedCycle,
out bool primeOrderFailed) || computedCycle == 0UL)
{
divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
divisor,
divisorData,
skipPrimeOrderHeuristic: primeOrderFailed);
}
else
{
divisorCycle = computedCycle;
}
if (divisorCycle == prime)
{
processedAll = true;
return true;
}
if (divisorCycle == 0UL)
{
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
}
}
remainingIterations--;
if (remainingIterations == 0UL)
{
processedAll = true;
return false;
}
divisor += step;
remainder10 = AddMod10(remainder10, step10);
remainder8 = AddMod8(remainder8, step8);
remainder3 = AddMod3(remainder3, step3);
remainder7 = AddMod7(remainder7, step7);
remainder11 = AddMod11(remainder11, step11);
}
}
private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
@@ -251,35 +390,67 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return residue == 1UL ? (byte)1 : (byte)0;
}
private static byte AddMod(byte value, byte delta, byte modulus)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod3(byte value, byte delta)
{
const int Modulus = 3;
int sum = value + delta;
if (sum >= modulus)
if (sum >= Modulus)
{
sum -= modulus;
sum -= Modulus;
}
return (byte)sum;
}
private void UpdateStatusUnsafe(ulong processedCount)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod7(byte value, byte delta)
{
if (processedCount == 0UL)
const int Modulus = 7;
int sum = value + delta;
if (sum >= Modulus)
{
return;
sum -= Modulus;
}
ulong interval = PerfectNumberConstants.ConsoleInterval;
if (interval == 0UL)
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)
{
_lastStatusDivisor = 0UL;
return;
sum -= Modulus;
}
ulong total = _lastStatusDivisor + processedCount;
// TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
// status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
_lastStatusDivisor = total % interval;
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 ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
-2
View File
@@ -227,7 +227,6 @@ public sealed class DivisorCycleCache
{
int length = divisors.Length;
var gpuLease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = gpuLease.EnterExecutionScope();
Accelerator accelerator = gpuLease.Accelerator;
var kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel);
@@ -309,7 +308,6 @@ public sealed class DivisorCycleCache
ArrayPool<byte>.Shared.Return(rentedStatus!, clearArray: false);
}
execution.Dispose();
gpuLease.Dispose();
}
+259 -20
View File
@@ -1,38 +1,277 @@
using System;
using System.Collections.Generic;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
internal static class DivisorGenerator
{
public static readonly uint[] SmallPrimes;
public static readonly ulong[] SmallPrimesPow2;
public static readonly uint[] SmallPrimesLastOne;
public static readonly ulong[] SmallPrimesPow2LastOne;
public static readonly uint[] SmallPrimesLastSeven;
public static readonly ulong[] SmallPrimesPow2LastSeven;
public static readonly uint[] SmallPrimesLastThree;
public static readonly ulong[] SmallPrimesPow2LastThree;
public static readonly uint[] SmallPrimesLastNine;
public static readonly ulong[] SmallPrimesPow2LastNine;
public static readonly uint[] SmallPrimesLastOneWithoutLastThree;
public static readonly ulong[] SmallPrimesPow2LastOneWithoutLastThree;
public static readonly uint[] SmallPrimesLastSevenWithoutLastThree;
public static readonly ulong[] SmallPrimesPow2LastSevenWithoutLastThree;
public static readonly uint[] SmallPrimesLastNineWithoutLastThree;
public static readonly ulong[] SmallPrimesPow2LastNineWithoutLastThree;
private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
private const ushort DecimalMaskWhenLastIsThree = (1 << 1) | (1 << 3) | (1 << 7) | (1 << 9);
private const ushort DecimalMaskWhenLastIsNine = (1 << 1) | (1 << 3) | (1 << 7) | (1 << 9);
private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsValidDivisor(
byte remainder10,
byte remainder8,
byte remainder3,
byte remainder5,
byte remainder7,
byte remainder11,
LastDigit lastDigit)
static DivisorGenerator()
{
ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
if (((decimalMask >> remainder10) & 1) == 0)
BuildSmallPrimes(
out SmallPrimes,
out SmallPrimesPow2,
out SmallPrimesLastOne,
out SmallPrimesPow2LastOne,
out SmallPrimesLastSeven,
out SmallPrimesPow2LastSeven,
out SmallPrimesLastThree,
out SmallPrimesPow2LastThree,
out SmallPrimesLastNine,
out SmallPrimesPow2LastNine,
out SmallPrimesLastOneWithoutLastThree,
out SmallPrimesPow2LastOneWithoutLastThree,
out SmallPrimesLastSevenWithoutLastThree,
out SmallPrimesPow2LastSevenWithoutLastThree,
out SmallPrimesLastNineWithoutLastThree,
out SmallPrimesPow2LastNineWithoutLastThree);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ushort GetDecimalMask(LastDigit lastDigit)
{
return lastDigit switch
{
return false;
LastDigit.Seven => DecimalMaskWhenLastIsSeven,
LastDigit.Three => DecimalMaskWhenLastIsThree,
LastDigit.Nine => DecimalMaskWhenLastIsNine,
_ => DecimalMaskOtherwise,
};
}
private static void BuildSmallPrimes(
out uint[] all,
out ulong[] allPow2,
out uint[] lastOne,
out ulong[] lastOnePow2,
out uint[] lastSeven,
out ulong[] lastSevenPow2,
out uint[] lastThree,
out ulong[] lastThreePow2,
out uint[] lastNine,
out ulong[] lastNinePow2,
out uint[] lastOneWithoutThree,
out ulong[] lastOneWithoutThreePow2,
out uint[] lastSevenWithoutThree,
out ulong[] lastSevenWithoutThreePow2,
out uint[] lastNineWithoutThree,
out ulong[] lastNineWithoutThreePow2)
{
uint target = PerfectNumberConstants.PrimesLimit;
int targetInt = checked((int)target);
all = new uint[targetInt];
allPow2 = new ulong[targetInt];
lastOne = new uint[targetInt];
lastOnePow2 = new ulong[targetInt];
lastSeven = new uint[targetInt];
lastSevenPow2 = new ulong[targetInt];
lastThree = new uint[targetInt];
lastThreePow2 = new ulong[targetInt];
lastNine = new uint[targetInt];
lastNinePow2 = new ulong[targetInt];
List<uint> primes = new(targetInt * 4 / 3 + 16);
List<uint> lastOneWithoutThreeList = new(targetInt);
List<ulong> lastOneWithoutThreePow2List = new(targetInt);
List<uint> lastSevenWithoutThreeList = new(targetInt);
List<ulong> lastSevenWithoutThreePow2List = new(targetInt);
List<uint> lastNineWithoutThreeList = new(targetInt);
List<ulong> lastNineWithoutThreePow2List = new(targetInt);
uint candidate = 2U;
int allCount = 0;
int lastOneCount = 0;
int lastSevenCount = 0;
int lastThreeCount = 0;
int lastNineCount = 0;
while (allCount < targetInt ||
lastOneCount < targetInt ||
lastSevenCount < targetInt ||
lastThreeCount < targetInt ||
lastNineCount < targetInt)
{
bool isPrime = true;
int primesCount = primes.Count;
for (int i = 0; i < primesCount; i++)
{
uint p = primes[i];
if (p * p > candidate)
{
break;
}
if (candidate % p == 0U)
{
isPrime = false;
break;
}
}
if (isPrime)
{
primes.Add(candidate);
ulong candidateSquared = checked((ulong)candidate * candidate);
if (allCount < targetInt)
{
all[allCount] = candidate;
allPow2[allCount] = candidateSquared;
allCount++;
}
if (candidate != 2U)
{
LastDigit lastDigit = GetLastDigit(candidate);
bool allowedForLastOne = IsAllowedForLastOne(candidate, lastDigit);
bool allowedForLastSeven = IsAllowedForLastSeven(candidate, lastDigit);
bool allowedForLastThree = IsAllowedForLastThree(candidate, lastDigit);
bool allowedForLastNine = IsAllowedForLastNine(candidate, lastDigit);
bool lastDigitIsThree = lastDigit == LastDigit.Three;
if (allowedForLastOne && lastOneCount < targetInt)
{
lastOne[lastOneCount] = candidate;
lastOnePow2[lastOneCount] = candidateSquared;
lastOneCount++;
}
if (allowedForLastSeven && lastSevenCount < targetInt)
{
lastSeven[lastSevenCount] = candidate;
lastSevenPow2[lastSevenCount] = candidateSquared;
lastSevenCount++;
}
if (allowedForLastThree && lastThreeCount < targetInt)
{
lastThree[lastThreeCount] = candidate;
lastThreePow2[lastThreeCount] = candidateSquared;
lastThreeCount++;
}
if (allowedForLastNine && lastNineCount < targetInt)
{
lastNine[lastNineCount] = candidate;
lastNinePow2[lastNineCount] = candidateSquared;
lastNineCount++;
}
if (!lastDigitIsThree && allowedForLastOne)
{
lastOneWithoutThreeList.Add(candidate);
lastOneWithoutThreePow2List.Add(candidateSquared);
}
if (!lastDigitIsThree && allowedForLastSeven)
{
lastSevenWithoutThreeList.Add(candidate);
lastSevenWithoutThreePow2List.Add(candidateSquared);
}
if (!lastDigitIsThree && allowedForLastNine)
{
lastNineWithoutThreeList.Add(candidate);
lastNineWithoutThreePow2List.Add(candidateSquared);
}
}
}
candidate = candidate == 2U ? 3U : candidate + 2U;
}
if (remainder8 != 1 && remainder8 != 7)
{
return false;
}
lastOneWithoutThree = lastOneWithoutThreeList.ToArray();
lastOneWithoutThreePow2 = lastOneWithoutThreePow2List.ToArray();
lastSevenWithoutThree = lastSevenWithoutThreeList.ToArray();
lastSevenWithoutThreePow2 = lastSevenWithoutThreePow2List.ToArray();
lastNineWithoutThree = lastNineWithoutThreeList.ToArray();
lastNineWithoutThreePow2 = lastNineWithoutThreePow2List.ToArray();
}
if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static LastDigit GetLastDigit(uint prime)
{
return (prime % 10U) switch
{
return false;
}
1U => LastDigit.One,
3U => LastDigit.Three,
5U => LastDigit.Five,
7U => LastDigit.Seven,
9U => LastDigit.Nine,
_ => throw new ArgumentOutOfRangeException(nameof(prime)),
};
}
return true;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastOne(uint prime, LastDigit lastDigit)
{
// TODO: Swap the `% 10` usage for ULongExtensions.Mod10 so the hot classification path
// reuses the benchmarked residue helper instead of repeated divisions.
return lastDigit switch
{
LastDigit.One or LastDigit.Three or LastDigit.Nine => true,
LastDigit.Seven => prime == 7U,
_ => prime == 11U,
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastSeven(uint prime, LastDigit lastDigit)
{
// TODO: Route this `% 10` classification through ULongExtensions.Mod10 to match the faster
// residue helper used elsewhere in the scanner.
return lastDigit switch
{
LastDigit.Three or LastDigit.Seven or LastDigit.Nine => true,
_ => prime == 11U,
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastThree(uint prime, LastDigit lastDigit)
{
// TODO: Replace the `% 10` residue classification with ULongExtensions.Mod10 once the helper
// migrates to UInt32 so we can reuse the optimized decimal residue path everywhere.
return lastDigit switch
{
LastDigit.Three or LastDigit.Seven => true,
LastDigit.One => prime == 11U,
LastDigit.Nine => prime == 19U,
_ => false,
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastNine(uint prime, LastDigit lastDigit)
{
// TODO: Replace the `% 10` residue classification with ULongExtensions.Mod10 when the helper
// supports UInt32 inputs so the residue tables share the optimized decimal path.
return lastDigit switch
{
LastDigit.Three or LastDigit.Seven or LastDigit.Nine => true,
LastDigit.One => prime == 11U,
_ => false,
};
}
}
@@ -74,7 +74,7 @@ internal struct ExponentRemainderStepperCpu
// single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the
// cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
_currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
PreviousExponent = exponent;
return ReduceCurrent();
}
@@ -92,7 +92,7 @@ internal struct ExponentRemainderStepperCpu
// TODO: Reuse the divisor-cycle derived Montgomery delta once the cache exposes single-cycle
// lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
_currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
PreviousExponent = exponent;
return _currentMontgomery == _montgomeryOne;
}
@@ -126,7 +126,7 @@ internal struct ExponentRemainderStepperCpu
// TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead
// of relying on the caller-provided delta so incremental scans remain in sync with the snapshot
// without computing additional cycles or mutating cache state.
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime);
_currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(montgomeryDelta, _modulus, _nPrime);
PreviousExponent = exponent;
isUnity = _currentMontgomery == _montgomeryOne;
return true;
@@ -138,5 +138,5 @@ internal struct ExponentRemainderStepperCpu
PreviousExponent = exponent;
}
private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);
private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime);
}
@@ -0,0 +1,142 @@
using System;
using System.Collections.Concurrent;
using System.Runtime.CompilerServices;
using ILGPU;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal static class DivisorGeneratorGpuCache
{
private sealed class AcceleratorReferenceComparer : IEqualityComparer<Accelerator>
{
public static AcceleratorReferenceComparer Instance { get; } = new();
public bool Equals(Accelerator? x, Accelerator? y)
{
return ReferenceEquals(x, y);
}
public int GetHashCode(Accelerator obj)
{
return RuntimeHelpers.GetHashCode(obj);
}
}
private sealed class DivisorGeneratorGpuTables : IDisposable
{
internal DivisorGeneratorGpuTables(Accelerator accelerator)
{
SmallPrimes = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimes);
SmallPrimesPow2 = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2);
SmallPrimesLastOne = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastOne);
SmallPrimesPow2LastOne = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastOne);
SmallPrimesLastSeven = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastSeven);
SmallPrimesPow2LastSeven = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastSeven);
SmallPrimesLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastThree);
SmallPrimesPow2LastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastThree);
SmallPrimesLastNine = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastNine);
SmallPrimesPow2LastNine = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastNine);
SmallPrimesLastOneWithoutLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastOneWithoutLastThree);
SmallPrimesPow2LastOneWithoutLastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree);
SmallPrimesLastSevenWithoutLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastSevenWithoutLastThree);
SmallPrimesPow2LastSevenWithoutLastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree);
SmallPrimesLastNineWithoutLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastNineWithoutLastThree);
SmallPrimesPow2LastNineWithoutLastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree);
}
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimes { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2 { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastOne { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastOne { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastSeven { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastSeven { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastThree { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastThree { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastNine { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastNine { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastOneWithoutLastThree { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastOneWithoutLastThree { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastSevenWithoutLastThree { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastSevenWithoutLastThree { get; }
internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastNineWithoutLastThree { get; }
internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastNineWithoutLastThree { get; }
public void Dispose()
{
SmallPrimes.Dispose();
SmallPrimesPow2.Dispose();
SmallPrimesLastOne.Dispose();
SmallPrimesPow2LastOne.Dispose();
SmallPrimesLastSeven.Dispose();
SmallPrimesPow2LastSeven.Dispose();
SmallPrimesLastThree.Dispose();
SmallPrimesPow2LastThree.Dispose();
SmallPrimesLastNine.Dispose();
SmallPrimesPow2LastNine.Dispose();
SmallPrimesLastOneWithoutLastThree.Dispose();
SmallPrimesPow2LastOneWithoutLastThree.Dispose();
SmallPrimesLastSevenWithoutLastThree.Dispose();
SmallPrimesPow2LastSevenWithoutLastThree.Dispose();
SmallPrimesLastNineWithoutLastThree.Dispose();
SmallPrimesPow2LastNineWithoutLastThree.Dispose();
}
private static MemoryBuffer1D<uint, Stride1D.Dense> AllocateAndCopyUInt(Accelerator accelerator, uint[] source)
{
var buffer = accelerator.Allocate1D<uint>(source.Length);
if (source.Length > 0)
{
buffer.View.CopyFromCPU(source);
}
return buffer;
}
private static MemoryBuffer1D<ulong, Stride1D.Dense> AllocateAndCopyULong(Accelerator accelerator, ulong[] source)
{
var buffer = accelerator.Allocate1D<ulong>(source.Length);
if (source.Length > 0)
{
buffer.View.CopyFromCPU(source);
}
return buffer;
}
}
private static readonly ConcurrentDictionary<Accelerator, DivisorGeneratorGpuTables> Cache = new(AcceleratorReferenceComparer.Instance);
public static void WarmUpByDivisorTables(Accelerator accelerator)
{
if (accelerator is null || accelerator.AcceleratorType == AcceleratorType.CPU)
{
return;
}
Cache.GetOrAdd(accelerator, acc => new DivisorGeneratorGpuTables(acc));
}
public static void Clear(Accelerator accelerator)
{
if (accelerator is null)
{
return;
}
if (Cache.TryRemove(accelerator, out var tables))
{
tables.Dispose();
}
}
public static void DisposeAll()
{
foreach (var entry in Cache)
{
entry.Value.Dispose();
}
Cache.Clear();
}
}
+52 -98
View File
@@ -1,3 +1,4 @@
using System;
using System.Collections.Concurrent;
using ILGPU;
using ILGPU.Runtime;
@@ -6,26 +7,19 @@ namespace PerfectNumbers.Core.Gpu;
public static class GpuContextPool
{
private static readonly bool PoolingEnabled = true;
// The pool intentionally skips pre-loading ProcessEightBitWindows kernels so accelerator initialization
// stays stable; GpuKernelPool JIT-compiles them on first use alongside the small-cycle uploads.
// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
public static bool ForceCpu { get; set; } = false;
// The pool intentionally skips pre-loading ProcessEightBitWindows kernels so accelerator initialization
// stays stable; GpuKernelPool JIT-compiles them on first use alongside the small-cycle uploads.
// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
internal sealed class PooledContext
{
public Context Context { get; }
public Accelerator Accelerator { get; }
internal sealed class PooledContext
{
private bool disposed;
public Context Context { get; }
public Accelerator Accelerator { get; }
public bool IsCpu { get; }
public object ExecutionLock { get; } = new();
public PooledContext(bool preferCpu)
public PooledContext()
{
Context = Context.CreateDefault();
Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context);
DivisorGeneratorGpuCache.WarmUpByDivisorTables(Accelerator);
// NOTE: Avoid loading/compiling any kernel here to prevent implicit
// CL stream/queue creation during accelerator construction.
// Some OpenCL drivers are fragile when a queue is created immediately
@@ -34,130 +28,90 @@ public static class GpuContextPool
public void Dispose()
{
if (disposed)
{
return;
}
// Ensure all cached GPU buffers for this accelerator are released.
NttGpuMath.ClearCaches(Accelerator);
// Release PrimeTester GPU state for this accelerator (kernel/device primes)
// Clear any per-accelerator cached resources to avoid releasing
// after the accelerator is destroyed.
DivisorGeneratorGpuCache.Clear(Accelerator);
PrimeTester.ClearGpuCaches(Accelerator);
Accelerator.Dispose();
Context.Dispose();
disposed = true;
}
}
private static readonly ConcurrentQueue<PooledContext> CpuPool = new();
private static readonly ConcurrentQueue<PooledContext> GpuPool = new();
private static readonly object CreationLock = new();
public static GpuContextLease Rent()
{
return RentPreferred(ForceCpu);
}
private static int WarmedGpuContextCount;
// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
public static GpuContextLease RentPreferred(bool preferCpu)
public static GpuContextLease Rent()
{
if (PoolingEnabled)
if (GpuPool.TryDequeue(out var gpu))
{
if (preferCpu && CpuPool.TryDequeue(out var cpu))
{
return new GpuContextLease(cpu, CreationLock);
}
if (!preferCpu && GpuPool.TryDequeue(out var gpu))
{
return new GpuContextLease(gpu, CreationLock);
}
return new GpuContextLease(gpu);
}
// Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
lock (CreationLock)
{
return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
}
// Create a new accelerator when the pool does not have one available.
return new GpuContextLease(new PooledContext());
}
public static void DisposeAll()
public static void WarmUpPool(int threadCount)
{
if (!PoolingEnabled)
if (threadCount <= WarmedGpuContextCount)
{
return;
}
while (CpuPool.TryDequeue(out var cpu))
int toCreate = threadCount - WarmedGpuContextCount;
var contexts = new PooledContext[toCreate];
for (int i = 0; i < toCreate; i++)
{
cpu.Dispose();
contexts[i] = new PooledContext();
}
for (int i = 0; i < toCreate; i++)
{
GpuPool.Enqueue(contexts[i]);
}
WarmedGpuContextCount = threadCount;
}
public static void DisposeAll()
{
while (GpuPool.TryDequeue(out var gpu))
{
gpu.Dispose();
}
DivisorGeneratorGpuCache.DisposeAll();
PrimeTester.DisposeGpuContexts();
WarmedGpuContextCount = 0;
}
private static void Return(PooledContext ctx)
{
if (PoolingEnabled)
ctx.Accelerator.Synchronize();
GpuPool.Enqueue(ctx);
}
public readonly struct GpuContextLease
{
private readonly PooledContext _ctx;
internal GpuContextLease(PooledContext ctx)
{
ctx.Accelerator.Synchronize();
if (ctx.IsCpu)
{
CpuPool.Enqueue(ctx);
}
else
{
GpuPool.Enqueue(ctx);
}
_ctx = ctx;
}
else
{
lock (CreationLock)
{
ctx.Dispose();
}
}
}
public struct GpuContextLease
{
private readonly PooledContext _ctx;
private bool disposed;
internal GpuContextLease(PooledContext ctx, object kernelInitLock)
{
_ctx = ctx;
KernelInitLock = kernelInitLock;
}
public Context Context => _ctx.Context;
public Accelerator Accelerator => _ctx.Accelerator;
public object KernelInitLock { get; }
public object ExecutionLock => _ctx.ExecutionLock;
public void Dispose()
{
if (disposed)
{
return;
}
public readonly Context Context => _ctx.Context;
public readonly Accelerator Accelerator => _ctx.Accelerator;
public readonly void Dispose()
{
Return(_ctx);
disposed = true;
}
}
}
}
}
+45 -105
View File
@@ -1,3 +1,4 @@
using System;
using System.Collections.Concurrent;
using ILGPU;
using ILGPU.Runtime;
@@ -9,30 +10,26 @@ public sealed class GpuKernelLease
{
private static readonly ConcurrentQueue<GpuKernelLease> Pool = new();
private IDisposable? _limiter;
private Action? _releaseLimiter;
private GpuContextLease? _gpu;
private KernelContainer? _kernels;
private AcceleratorStream? _stream;
private object? _executionLock;
private bool _disposed;
private GpuKernelLease()
{
}
internal static GpuKernelLease Rent(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels)
internal static GpuKernelLease Rent(GpuPrimeWorkLimiter.Lease limiter, GpuContextLease gpu, KernelContainer kernels)
{
if (!Pool.TryDequeue(out var lease))
{
lease = new GpuKernelLease();
}
lease._limiter = limiter;
lease._releaseLimiter = limiter.Dispose;
lease._gpu = gpu;
lease._kernels = kernels;
lease._stream = null;
lease._executionLock = gpu.ExecutionLock;
lease._disposed = false;
return lease;
}
@@ -54,22 +51,17 @@ public sealed class GpuKernelLease
}
}
public ExecutionScope EnterExecutionScope() => new ExecutionScope(_executionLock);
public Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>> OrderKernel
{
get
{
var accel = Accelerator; // avoid capturing 'this' in lambda
lock (_gpu!.Value.KernelInitLock)
return KernelContainer.InitOnce(ref _kernels!.Order, () =>
{
return KernelContainer.InitOnce(ref _kernels!.Order, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernels.OrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>>();
});
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernels.OrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>>();
});
}
}
@@ -78,15 +70,12 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
return KernelContainer.InitOnce(ref _kernels!.Pow2Mod, () =>
{
return KernelContainer.InitOnce(ref _kernels!.Pow2Mod, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernels.Pow2ModKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernels.Pow2ModKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
@@ -95,15 +84,12 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
return KernelContainer.InitOnce(ref _kernels!.Incremental, () =>
{
return KernelContainer.InitOnce(ref _kernels!.Incremental, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernels.IncrementalKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernels.IncrementalKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
@@ -112,15 +98,12 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
return KernelContainer.InitOnce(ref _kernels!.IncrementalOrder, () =>
{
return KernelContainer.InitOnce(ref _kernels!.IncrementalOrder, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernels.IncrementalOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernels.IncrementalOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
@@ -129,15 +112,12 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
return KernelContainer.InitOnce(ref _kernels!.Pow2ModOrder, () =>
{
return KernelContainer.InitOnce(ref _kernels!.Pow2ModOrder, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernels.Pow2ModOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernels.Pow2ModOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
@@ -146,68 +126,28 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
return KernelContainer.InitOnce(ref _kernels!.SmallPrimeFactor, () =>
{
return KernelContainer.InitOnce(ref _kernels!.SmallPrimeFactor, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
private void Dispose(bool disposing)
{
if (_disposed)
{
return;
}
if (disposing)
{
_stream?.Dispose();
_stream = null;
_gpu?.Dispose();
_gpu = null;
_executionLock = null;
_limiter?.Dispose();
_limiter = null;
}
_kernels = null;
_disposed = true;
Pool.Enqueue(this);
}
public void Dispose()
{
Dispose(disposing: false);
_stream?.Dispose();
_stream = null;
_gpu?.Dispose();
_gpu = null;
_releaseLimiter?.Invoke();
_releaseLimiter = null;
_kernels = null;
Pool.Enqueue(this);
}
public readonly struct ExecutionScope
{
private readonly object? _lock;
public ExecutionScope(object? sync)
{
_lock = sync;
if (_lock is not null)
{
Monitor.Enter(_lock);
}
}
public void Dispose()
{
if (_lock is not null)
{
Monitor.Exit(_lock);
}
}
}
}
+1 -3
View File
@@ -280,7 +280,7 @@ public class GpuKernelPool
public static GpuKernelLease GetKernel(bool useGpuOrder)
{
var limiter = GpuPrimeWorkLimiter.Acquire();
var gpu = RentPreferred(preferCpu: !useGpuOrder);
var gpu = Rent();
var accelerator = gpu.Accelerator;
var kernels = GetKernels(accelerator);
return GpuKernelLease.Rent(limiter, gpu, kernels);
@@ -293,12 +293,10 @@ public class GpuKernelPool
public static void Run(Action<Accelerator, AcceleratorStream> action)
{
var lease = GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
var accelerator = lease.Accelerator;
var stream = accelerator.CreateStream();
action(accelerator, stream);
stream.Dispose();
execution.Dispose();
lease.Dispose();
}
+15 -31
View File
@@ -3,58 +3,42 @@ namespace PerfectNumbers.Core.Gpu;
public static class GpuPrimeWorkLimiter
{
private static SemaphoreSlim _semaphore = new(1, int.MaxValue);
private static readonly object _lock = new();
private static int _currentLimit = 1;
public static IDisposable Acquire()
public static Lease Acquire()
{
// TODO: Inline the pooled limiter guard from the GpuLimiterThroughputBenchmarks so prime-sieve
// GPU jobs stop allocating Releaser instances on every acquisition.
// GPU jobs stop allocating Lease instances on every acquisition.
var sem = _semaphore;
sem.Wait();
return new Releaser(sem);
// sem.Wait();
return new Lease(sem);
}
public static void SetLimit(int value)
{
if (value < 1)
{
value = 1;
}
if (value == _currentLimit)
{
return;
}
lock (_lock)
{
if (value == _currentLimit)
{
return;
}
// TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using
// the same semaphore pool without reallocating per adjustment.
_semaphore = new SemaphoreSlim(value, value);
_currentLimit = value;
}
// TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using
// the same semaphore pool without reallocating per adjustment.
_semaphore = new SemaphoreSlim(value, value);
_currentLimit = value;
}
private sealed class Releaser : IDisposable
public sealed class Lease
{
private readonly SemaphoreSlim _sem;
private bool _disposed;
public Releaser(SemaphoreSlim sem)
internal Lease(SemaphoreSlim sem)
{
_sem = sem;
}
public void Dispose()
{
if (_disposed)
{
return;
}
_sem.Release();
_disposed = true;
// _sem.Release();
}
}
}
+6 -13
View File
@@ -6,14 +6,14 @@ public static class GpuWorkLimiter
private static readonly object _lock = new();
private static int _currentLimit = 1;
public static IDisposable Acquire()
public static Lease Acquire()
{
// TODO: Replace the per-call Releaser allocation with the pooled struct-based guard from the
// TODO: Replace the per-call Lease allocation with the pooled struct-based guard from the
// GpuLimiterThroughputBenchmarks so limiter acquisition does not allocate when we enter the
// GPU scanning hot path.
var sem = _semaphore;
sem.Wait();
return new Releaser(sem);
return new Lease(sem);
}
public static void SetLimit(int value)
@@ -32,7 +32,7 @@ public static class GpuWorkLimiter
// Swap to a new semaphore for the new limit. We intentionally
// do not dispose the old semaphore to avoid releasing handles
// that may still be in use by existing Releaser instances.
// that may still be in use by existing Lease instances.
// TODO: Consolidate with GpuPrimeWorkLimiter so both limiters share a pooled SemaphoreSlim and avoid rebuilding the
// limiter state whenever limits are adjusted from the CLI.
_semaphore = new SemaphoreSlim(value, value);
@@ -40,25 +40,18 @@ public static class GpuWorkLimiter
}
}
private sealed class Releaser : IDisposable
public sealed class Lease
{
private readonly SemaphoreSlim _sem;
private bool _disposed;
public Releaser(SemaphoreSlim sem)
internal Lease(SemaphoreSlim sem)
{
_sem = sem;
}
public void Dispose()
{
if (_disposed)
{
return;
}
_sem.Release();
_disposed = true;
}
}
}
@@ -0,0 +1,97 @@
using System.Runtime.CompilerServices;
using ILGPU;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal readonly struct HeuristicGpuDivisorTables
{
public static readonly HeuristicGpuDivisorTables Shared = default;
public readonly ArrayView1D<ulong, Stride1D.Dense> GroupADivisors;
public readonly ArrayView1D<ulong, Stride1D.Dense> GroupBDivisorsEnding1;
public readonly ArrayView1D<ulong, Stride1D.Dense> GroupBDivisorsEnding7;
public readonly ArrayView1D<ulong, Stride1D.Dense> GroupBDivisorsEnding9;
public HeuristicGpuDivisorTables(
ArrayView1D<ulong, Stride1D.Dense> groupADivisors,
ArrayView1D<ulong, Stride1D.Dense> groupBDivisorsEnding1,
ArrayView1D<ulong, Stride1D.Dense> groupBDivisorsEnding7,
ArrayView1D<ulong, Stride1D.Dense> groupBDivisorsEnding9)
{
GroupADivisors = groupADivisors;
GroupBDivisorsEnding1 = groupBDivisorsEnding1;
GroupBDivisorsEnding7 = groupBDivisorsEnding7;
GroupBDivisorsEnding9 = groupBDivisorsEnding9;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public ArrayView1D<ulong, Stride1D.Dense> SelectGroupB(byte ending)
{
return ending switch
{
1 => GroupBDivisorsEnding1,
7 => GroupBDivisorsEnding7,
9 => GroupBDivisorsEnding9,
_ => ArrayView1D<ulong, Stride1D.Dense>.Empty,
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public ArrayView1D<ulong, Stride1D.Dense> SelectDivisors(HeuristicGpuDivisorTableKind kind, byte ending)
{
return kind switch
{
HeuristicGpuDivisorTableKind.GroupA => GroupADivisors,
HeuristicGpuDivisorTableKind.GroupBEnding1 => GroupBDivisorsEnding1,
HeuristicGpuDivisorTableKind.GroupBEnding7 => GroupBDivisorsEnding7,
HeuristicGpuDivisorTableKind.GroupBEnding9 => GroupBDivisorsEnding9,
_ => ArrayView1D<ulong, Stride1D.Dense>.Empty,
};
}
internal static void InitializeShared(in HeuristicGpuDivisorTables tables)
{
Unsafe.AsRef(in Shared) = tables;
}
}
internal readonly struct HeuristicGpuCombinedDivisorTables
{
public static readonly HeuristicGpuCombinedDivisorTables Shared = default;
public readonly ArrayView1D<ulong, Stride1D.Dense> CombinedDivisorsEnding1;
public readonly ArrayView1D<ulong, Stride1D.Dense> CombinedDivisorsEnding3;
public readonly ArrayView1D<ulong, Stride1D.Dense> CombinedDivisorsEnding7;
public readonly ArrayView1D<ulong, Stride1D.Dense> CombinedDivisorsEnding9;
public HeuristicGpuCombinedDivisorTables(
ArrayView1D<ulong, Stride1D.Dense> combinedDivisorsEnding1,
ArrayView1D<ulong, Stride1D.Dense> combinedDivisorsEnding3,
ArrayView1D<ulong, Stride1D.Dense> combinedDivisorsEnding7,
ArrayView1D<ulong, Stride1D.Dense> combinedDivisorsEnding9)
{
CombinedDivisorsEnding1 = combinedDivisorsEnding1;
CombinedDivisorsEnding3 = combinedDivisorsEnding3;
CombinedDivisorsEnding7 = combinedDivisorsEnding7;
CombinedDivisorsEnding9 = combinedDivisorsEnding9;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public ArrayView1D<ulong, Stride1D.Dense> SelectCombined(byte ending)
{
return ending switch
{
1 => CombinedDivisorsEnding1,
3 => CombinedDivisorsEnding3,
7 => CombinedDivisorsEnding7,
9 => CombinedDivisorsEnding9,
_ => ArrayView1D<ulong, Stride1D.Dense>.Empty,
};
}
internal static void InitializeShared(in HeuristicGpuCombinedDivisorTables tables)
{
Unsafe.AsRef(in Shared) = tables;
}
}
@@ -36,30 +36,45 @@ internal static class DivisorByDivisorKernels
byte firstHitValue;
ulong previousExponent = firstExponent;
ulong cycleRemainder = 0UL;
if (cycleLength == 0UL)
{
firstHitValue = firstUnity ? (byte)1 : (byte)0;
}
else
bool hasCycle = cycleLength != 0UL;
if (hasCycle)
{
cycleRemainder = firstExponent % cycleLength;
firstHitValue = cycleRemainder == 0UL && firstUnity ? (byte)1 : (byte)0;
}
else
{
firstHitValue = firstUnity ? (byte)1 : (byte)0;
}
hits[offset] = firstHitValue;
byte anyHit = firstHitValue;
for (int i = 1; i < count; i++)
int endIndex = offset + count;
if (hasCycle)
{
ulong exponent = exponents[offset + i];
ulong delta = exponent - previousExponent;
previousExponent = exponent;
if (cycleLength != 0UL)
for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++)
{
var remainderValue = new GpuUInt128(cycleRemainder);
remainderValue.AddMod(delta, cycleLength);
cycleRemainder = remainderValue.Low;
ulong exponent = exponents[exponentIndex];
ulong delta = exponent - previousExponent;
previousExponent = exponent;
if (delta >= cycleLength)
{
delta -= cycleLength;
if (delta >= cycleLength)
{
delta %= cycleLength;
}
}
ulong nextRemainder = cycleRemainder + delta;
if (nextRemainder >= cycleLength)
{
nextRemainder -= cycleLength;
}
cycleRemainder = nextRemainder;
byte hit = 0;
if (cycleRemainder == 0UL)
@@ -67,13 +82,17 @@ internal static class DivisorByDivisorKernels
hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0;
}
hits[offset + i] = hit;
hits[exponentIndex] = hit;
anyHit |= hit;
}
else
}
else
{
for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++)
{
ulong exponent = exponents[exponentIndex];
byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0;
hits[offset + i] = hit;
hits[exponentIndex] = hit;
anyHit |= hit;
}
}
@@ -1,33 +1,73 @@
using System.Runtime.CompilerServices;
using ILGPU;
using ILGPU.Algorithms;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal enum HeuristicGpuDivisorTableKind : byte
{
GroupA = 0,
GroupBEnding1 = 1,
GroupBEnding7 = 7,
GroupBEnding9 = 9,
Combined = 255,
}
internal static class PrimeTesterKernels
{
// GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise.
public static void SmallPrimeSieveKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<uint> smallPrimes, ArrayView<byte> results)
public static void SmallPrimeSieveKernel(
Index1D index,
ArrayView<ulong> numbers,
ArrayView<uint> smallPrimes,
ArrayView<uint> smallPrimesLastOne,
ArrayView<uint> smallPrimesLastSeven,
ArrayView<uint> smallPrimesLastThree,
ArrayView<uint> smallPrimesLastNine,
ArrayView<ulong> smallPrimesPow2,
ArrayView<ulong> smallPrimesPow2LastOne,
ArrayView<ulong> smallPrimesPow2LastSeven,
ArrayView<ulong> smallPrimesPow2LastThree,
ArrayView<ulong> smallPrimesPow2LastNine,
ArrayView<byte> results)
{
ulong n = numbers[index];
byte result = 1;
// EvenPerfectBitScanner filters candidates so production GPU launches only see:
// - n >= 31 (the scanner starts at 136,279,841 and test mode seeds 31)
// - n odd (candidate generators and ModResidueTracker eliminate even exponents)
// - n not divisible by 5 (residue prefiltering and AddPrimes transform cover these composites)
// Tests should honor the same constraints so the host wrapper can remain branchless.
ArrayView<uint> primes = smallPrimes;
ArrayView<ulong> primeSquares = smallPrimesPow2;
ulong lastDigit = n % 10UL;
switch (lastDigit)
{
case 1UL:
primes = smallPrimesLastOne;
primeSquares = smallPrimesPow2LastOne;
break;
case 3UL:
primes = smallPrimesLastThree;
primeSquares = smallPrimesPow2LastThree;
break;
case 7UL:
primes = smallPrimesLastSeven;
primeSquares = smallPrimesPow2LastSeven;
break;
case 9UL:
primes = smallPrimesLastNine;
primeSquares = smallPrimesPow2LastNine;
break;
}
long length = smallPrimes.Length;
int length = (int)primes.Length;
for (int i = 0; i < length; i++)
{
ulong prime = smallPrimes[i];
ulong primeSquare = (ulong)prime * prime;
ulong primeSquare = primeSquares[i];
if (primeSquare > n)
{
break;
}
ulong prime = primes[i];
if (n % prime == 0UL)
{
result = 0;
@@ -38,16 +78,32 @@ internal static class PrimeTesterKernels
results[index] = result;
}
public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView<ulong> divisors, ulong n, ArrayView<byte> results)
public static void HeuristicTrialDivisionKernel(
Index1D index,
ArrayView<int> resultFlag,
ulong n,
ulong maxDivisor,
HeuristicGpuDivisorTableKind tableKind,
HeuristicGpuDivisorTables tables,
HeuristicGpuCombinedDivisorTables combinedTables)
{
ulong divisor = divisors[index];
if (divisor <= 1UL)
{
results[index] = 0;
return;
}
byte nMod10 = (byte)(n % 10UL);
ArrayView1D<ulong, Stride1D.Dense> divisors = tableKind == HeuristicGpuDivisorTableKind.Combined
? combinedTables.SelectCombined(nMod10)
: tables.SelectDivisors(tableKind, nMod10);
results[index] = n % divisor == 0UL ? (byte)1 : (byte)0;
int divisorLength = (int)divisors.Length;
int threadIndex = index;
bool skipThread = threadIndex >= divisorLength;
ulong divisor = skipThread ? 1UL : divisors[threadIndex];
skipThread = skipThread | (divisor > maxDivisor);
ulong remainder = skipThread ? 1UL : n % divisor;
bool hasFactor = !skipThread & (remainder == 0UL);
if (hasFactor)
{
resultFlag[0] = 1;
}
}
public static void SharesFactorKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<byte> results)
@@ -94,7 +94,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
var gpuLease = RentAccelerator();
var accelerator = gpuLease.Accelerator;
Monitor.Enter(gpuLease.ExecutionLock);
// Monitor.Enter(gpuLease.ExecutionLock);
var kernel = GetKernel(accelerator);
var resources = RentBatchResources(accelerator, batchCapacity);
@@ -112,7 +112,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
resources.HitIndexBuffer,
resources.Divisors,
resources.Exponents,
resources.FilteredDivisors,
resources.DivisorData,
resources.Offsets,
resources.Counts,
@@ -122,7 +121,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
out processedCount);
ReturnBatchResources(resources);
Monitor.Exit(gpuLease.ExecutionLock);
// Monitor.Exit(gpuLease.ExecutionLock);
ReturnAccelerator(gpuLease);
if (composite)
@@ -177,7 +176,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
MemoryBuffer1D<int, Stride1D.Dense> hitIndexBuffer,
ulong[] divisors,
ulong[] exponents,
ulong[] filteredDivisors,
GpuDivisorPartialData[] divisorData,
int[] offsets,
int[] counts,
@@ -192,7 +190,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
processedCount = 0UL;
lastProcessed = 0UL;
int chunkCapacity = Math.Max(1, batchCapacity);
// EvenPerfectBitScanner always configures non-empty GPU divisor buffers, so the defensive guard stays
// commented out. Tests and benchmarks that rely on zero-capacity buffers must be updated instead of re-enabling
// this branch.
// if (batchCapacity <= 0)
// {
// throw new InvalidOperationException("Divisor buffers must be non-empty.");
// }
int chunkCountBaseline = batchCapacity;
UInt128 twoP128 = (UInt128)prime << 1;
UInt128 allowedMax128 = allowedMax;
@@ -209,77 +215,103 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
ulong currentK = 1UL;
ulong remainingCount = maxK;
byte step10 = (byte)(twoP128 % 10UL);
byte step8 = (byte)(twoP128 % 8UL);
byte step5 = (byte)(twoP128 % 5UL);
byte step3 = (byte)(twoP128 % 3UL);
byte step7 = (byte)(twoP128 % 7UL);
byte step11 = (byte)(twoP128 % 11UL);
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);
UInt128 currentDivisor128 = firstDivisor128;
byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL));
byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL));
byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL));
byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL));
byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL));
byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL));
ulong currentDivisorValue = (ulong)currentDivisor128;
byte remainder10 = (byte)(currentDivisorValue % 10UL);
byte remainder8 = (byte)(currentDivisorValue & 7UL);
byte remainder3 = (byte)(currentDivisorValue % 3UL);
byte remainder7 = (byte)(currentDivisorValue % 7UL);
byte remainder11 = (byte)(currentDivisorValue % 11UL);
LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);
Span<ulong> filteredStorage = filteredDivisors.AsSpan();
Span<ulong> divisorStorage = divisors.AsSpan();
Span<ulong> exponentStorage = exponents.AsSpan();
Span<GpuDivisorPartialData> divisorDataStorage = divisorData.AsSpan();
Span<int> offsetStorage = offsets.AsSpan();
Span<int> countStorage = counts.AsSpan();
Span<ulong> cycleStorage = cycles.AsSpan();
while (currentK <= maxK && !composite)
Span<ulong> divisorSpan = divisorStorage;
Span<ulong> exponentSpan = exponentStorage;
Span<GpuDivisorPartialData> divisorDataSpan = divisorDataStorage;
Span<int> offsetSpan = offsetStorage;
Span<int> countSpan = countStorage;
Span<ulong> cycleSpan = cycleStorage;
ref GpuDivisorPartialData divisorDataRef = ref MemoryMarshal.GetReference(divisorDataSpan);
ref int offsetRef = ref MemoryMarshal.GetReference(offsetSpan);
ref int countRef = ref MemoryMarshal.GetReference(countSpan);
ref ulong cycleRef = ref MemoryMarshal.GetReference(cycleSpan);
ref ulong exponentRef = ref MemoryMarshal.GetReference(exponentSpan);
ArrayView1D<GpuDivisorPartialData, Stride1D.Dense> divisorDataView = divisorDataBuffer.View;
ArrayView1D<int, Stride1D.Dense> offsetView = offsetBuffer.View;
ArrayView1D<int, Stride1D.Dense> countView = countBuffer.View;
ArrayView1D<ulong, Stride1D.Dense> cycleView = cycleBuffer.View;
ArrayView1D<ulong, Stride1D.Dense> exponentViewDevice = exponentBuffer.View;
ArrayView1D<byte, Stride1D.Dense> hitsView = hitsBuffer.View;
ArrayView1D<int, Stride1D.Dense> hitIndexView = hitIndexBuffer.View;
while (remainingCount > 0UL && !composite)
{
int chunkCount = Math.Min(chunkCapacity, batchCapacity);
ulong remainingK = maxK - currentK + 1UL;
if ((ulong)chunkCount > remainingK)
int chunkCount = chunkCountBaseline;
if ((ulong)chunkCount > remainingCount)
{
chunkCount = (int)remainingK;
chunkCount = (int)remainingCount;
}
if (chunkCount <= 0)
{
processedAll = true;
break;
}
int filteredCount = 0;
UInt128 nextDivisor128 = currentDivisor128;
int admissibleCount = 0;
byte localRemainder10 = remainder10;
byte localRemainder8 = remainder8;
byte localRemainder5 = remainder5;
byte localRemainder3 = remainder3;
byte localRemainder7 = remainder7;
byte localRemainder11 = remainder11;
for (int i = 0; i < chunkCount; i++)
{
ulong candidate = (ulong)nextDivisor128;
if (DivisorGenerator.IsValidDivisor(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit))
bool passesSmallModuli = localRemainder3 != 0 && localRemainder7 != 0 && localRemainder11 != 0;
if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0)
{
filteredStorage[filteredCount++] = candidate;
ulong divisorValue = (ulong)nextDivisor128;
MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData);
if (divisorCycle == prime)
{
processedCount += (ulong)(i + 1);
lastProcessed = divisorValue;
coveredRange = true;
return true;
}
divisorSpan[admissibleCount] = divisorValue;
divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue);
offsetSpan[admissibleCount] = admissibleCount;
countSpan[admissibleCount] = 1;
cycleSpan[admissibleCount] = divisorCycle;
exponentSpan[admissibleCount] = prime;
admissibleCount++;
}
nextDivisor128 += twoP128;
localRemainder10 = AddMod(localRemainder10, step10, 10);
localRemainder8 = AddMod(localRemainder8, step8, 8);
localRemainder5 = AddMod(localRemainder5, step5, 5);
localRemainder3 = AddMod(localRemainder3, step3, 3);
localRemainder7 = AddMod(localRemainder7, step7, 7);
localRemainder11 = AddMod(localRemainder11, step11, 11);
localRemainder10 = AddMod10(localRemainder10, step10);
localRemainder8 = AddMod8(localRemainder8, step8);
localRemainder3 = AddMod3(localRemainder3, step3);
localRemainder7 = AddMod7(localRemainder7, step7);
localRemainder11 = AddMod11(localRemainder11, step11);
}
remainder10 = localRemainder10;
remainder8 = localRemainder8;
remainder5 = localRemainder5;
remainder3 = localRemainder3;
remainder7 = localRemainder7;
remainder11 = localRemainder11;
@@ -289,69 +321,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
UInt128 lastDivisor128 = nextDivisor128 - twoP128;
lastProcessed = (ulong)lastDivisor128;
currentDivisor128 = nextDivisor128;
currentK += (ulong)chunkCount;
processedAll = currentK > maxK;
if (filteredCount == 0)
{
continue;
}
Span<ulong> filteredDivisorsSpan = filteredStorage[..filteredCount];
Span<ulong> divisorSpan = divisorStorage;
Span<ulong> exponentSpan = exponentStorage;
Span<GpuDivisorPartialData> divisorDataSpan = divisorDataStorage;
Span<int> offsetSpan = offsetStorage;
Span<int> countSpan = countStorage;
Span<ulong> cycleSpan = cycleStorage;
int admissibleCount = 0;
int exponentIndex = 0;
for (int i = 0; i < filteredCount; i++)
{
ulong divisorValue = filteredDivisorsSpan[i];
MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData);
if (divisorCycle != prime)
{
continue;
}
divisorSpan[admissibleCount] = divisorValue;
exponentSpan[exponentIndex] = prime;
divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue);
offsetSpan[admissibleCount] = exponentIndex;
countSpan[admissibleCount] = 1;
cycleSpan[admissibleCount] = divisorCycle;
admissibleCount++;
exponentIndex++;
}
remainingCount -= (ulong)chunkCount;
if (admissibleCount == 0)
{
continue;
}
ArrayView1D<GpuDivisorPartialData, Stride1D.Dense> divisorView = divisorDataBuffer.View.SubView(0, admissibleCount);
ArrayView1D<int, Stride1D.Dense> offsetView = offsetBuffer.View.SubView(0, admissibleCount);
ArrayView1D<int, Stride1D.Dense> countView = countBuffer.View.SubView(0, admissibleCount);
ArrayView1D<ulong, Stride1D.Dense> cycleView = cycleBuffer.View.SubView(0, admissibleCount);
ArrayView1D<ulong, Stride1D.Dense> exponentView = exponentBuffer.View.SubView(0, exponentIndex);
ArrayView1D<byte, Stride1D.Dense> hitsView = hitsBuffer.View.SubView(0, exponentIndex);
ArrayView1D<int, Stride1D.Dense> hitIndexView = hitIndexBuffer.View.SubView(0, 1);
divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), admissibleCount);
offsetView.CopyFromCPU(ref MemoryMarshal.GetReference(offsetSpan), admissibleCount);
countView.CopyFromCPU(ref MemoryMarshal.GetReference(countSpan), admissibleCount);
cycleView.CopyFromCPU(ref MemoryMarshal.GetReference(cycleSpan), admissibleCount);
exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), exponentIndex);
divisorDataView.CopyFromCPU(ref divisorDataRef, admissibleCount);
offsetView.CopyFromCPU(ref offsetRef, admissibleCount);
countView.CopyFromCPU(ref countRef, admissibleCount);
cycleView.CopyFromCPU(ref cycleRef, admissibleCount);
exponentViewDevice.CopyFromCPU(ref exponentRef, admissibleCount);
int sentinel = int.MaxValue;
hitIndexView.CopyFromCPU(ref sentinel, 1);
kernel(admissibleCount, divisorView, offsetView, countView, exponentView, cycleView, hitsView, hitIndexView);
kernel(admissibleCount, divisorDataView, offsetView, countView, exponentViewDevice, cycleView, hitsView, hitIndexView);
int firstHit = sentinel;
hitIndexView.CopyToCPU(ref firstHit, 1);
@@ -362,6 +348,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex];
}
processedAll = remainingCount == 0UL;
coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128);
return composite;
}
@@ -379,23 +366,67 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod(byte value, byte delta, byte modulus)
private static byte AddMod3(byte value, byte delta)
{
int result = value + delta;
if (result >= modulus)
const int Modulus = 3;
int sum = value + delta;
if (sum >= Modulus)
{
result -= modulus;
if (result >= modulus)
{
result -= modulus;
}
sum -= Modulus;
}
return (byte)result;
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)
@@ -421,13 +452,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
// {
// return 0UL;
// }
ulong cappedLimit = divisorLimit;
const ulong TestPrimeDivisorCap = 170_000_000UL;
if (prime < 10_000_000UL && cappedLimit > TestPrimeDivisorCap)
{
cappedLimit = TestPrimeDivisorCap;
}
if (prime - 1UL >= 64UL)
{
return divisorLimit;
return cappedLimit;
}
ulong computedLimit = (1UL << (int)(prime - 1UL)) - 1UL;
return computedLimit < divisorLimit ? computedLimit : divisorLimit;
if (prime < 10_000_000UL && computedLimit > cappedLimit)
{
computedLimit = cappedLimit;
}
return computedLimit < cappedLimit ? computedLimit : cappedLimit;
}
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -445,19 +488,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private MemoryBuffer1D<byte, Stride1D.Dense> _hitBuffer = null!;
private ulong[] _hostBuffer = null!;
private int _capacity;
private bool _disposed;
internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner)
{
_owner = owner;
_lease = GpuContextPool.RentPreferred(preferCpu: false);
_lease = GpuContextPool.Rent();
_accelerator = _lease.Accelerator;
_capacity = Math.Max(1, owner._gpuBatchSize);
}
internal void Reset()
{
_disposed = false;
}
private void EnsureExecutionResourcesLocked(int requiredCapacity)
@@ -504,11 +545,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan<ulong> primes, Span<byte> hits)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(DivisorScanSession));
}
int length = primes.Length;
// EvenPerfectBitScanner always supplies at least one exponent per divisor check, so the guard stays commented out.
// if (length == 0)
@@ -544,7 +580,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
}
Monitor.Enter(_lease.ExecutionLock);
// Monitor.Enter(_lease.ExecutionLock);
EnsureExecutionResourcesLocked(length);
@@ -557,51 +593,39 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ArrayView1D<byte, Stride1D.Dense> hitView = _hitBuffer.View;
GpuDivisorPartialData partialData = new GpuDivisorPartialData(divisor);
divisorView.SubView(0, 1).CopyFromCPU(ref partialData, 1);
divisorView.CopyFromCPU(ref partialData, 1);
int offsetValue = 0;
offsetView.SubView(0, 1).CopyFromCPU(ref offsetValue, 1);
offsetView.CopyFromCPU(ref offsetValue, 1);
int countValue = length;
countView.SubView(0, 1).CopyFromCPU(ref countValue, 1);
countView.CopyFromCPU(ref countValue, 1);
ulong cycleValue = cycle;
cycleView.SubView(0, 1).CopyFromCPU(ref cycleValue, 1);
cycleView.CopyFromCPU(ref cycleValue, 1);
Span<ulong> hostSpan = _hostBuffer.AsSpan(0, length);
primes.CopyTo(hostSpan);
ref ulong hostRef = ref MemoryMarshal.GetReference(hostSpan);
ArrayView1D<ulong, Stride1D.Dense> exponentSlice = exponentView.SubView(0, length);
exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length);
ArrayView1D<byte, Stride1D.Dense> hitSliceView = hitView.SubView(0, length);
exponentView.CopyFromCPU(ref hostRef, length);
int sentinel = int.MaxValue;
firstHitView.SubView(0, 1).CopyFromCPU(ref sentinel, 1);
firstHitView.CopyFromCPU(ref sentinel, 1);
var kernel = _kernel;
if (kernel is null)
{
Monitor.Exit(_lease.ExecutionLock);
throw new InvalidOperationException("GPU divisor kernel was not initialized.");
}
var kernel = _kernel!;
kernel(1, divisorView, offsetView, countView, exponentSlice, cycleView, hitSliceView, firstHitView);
kernel(1, divisorView, offsetView, countView, exponentView, cycleView, hitView, firstHitView);
Span<byte> hitSlice = hits.Slice(0, length);
hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), length);
ref byte hitRef = ref MemoryMarshal.GetReference(hitSlice);
hitView.CopyToCPU(ref hitRef, length);
Monitor.Exit(_lease.ExecutionLock);
// Monitor.Exit(_lease.ExecutionLock);
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_owner._sessionPool.Add(this);
}
}
@@ -633,7 +657,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
return lease;
}
return GpuContextPool.RentPreferred(preferCpu: false);
return GpuContextPool.Rent();
}
private void ReturnAccelerator(GpuContextPool.GpuContextLease lease) => _acceleratorPool.Add(lease);
@@ -651,7 +675,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private void ReturnBatchResources(BatchResources resources) => _resourcePool.Add(resources);
private sealed class BatchResources : IDisposable
private sealed class BatchResources
{
private Accelerator? _accelerator;
private MemoryBuffer1D<GpuDivisorPartialData, Stride1D.Dense>? _divisorDataBuffer;
@@ -688,8 +712,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
internal ulong[] Exponents { get; private set; } = null!;
internal ulong[] FilteredDivisors { get; private set; } = null!;
internal GpuDivisorPartialData[] DivisorData { get; private set; } = null!;
internal int[] Offsets { get; private set; } = null!;
@@ -764,7 +786,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
Divisors = ulongPool.Rent(capacity);
Exponents = ulongPool.Rent(capacity);
FilteredDivisors = ulongPool.Rent(capacity);
DivisorData = partialPool.Rent(capacity);
Offsets = intPool.Rent(capacity);
Counts = intPool.Rent(capacity);
@@ -803,14 +824,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ArrayPool<GpuDivisorPartialData> partialPool = ThreadStaticPools.GpuDivisorPartialDataPool;
ulongPool.Return(Divisors, clearArray: false);
ulongPool.Return(Exponents, clearArray: false);
ulongPool.Return(FilteredDivisors, clearArray: false);
partialPool.Return(DivisorData, clearArray: false);
intPool.Return(Offsets, clearArray: false);
intPool.Return(Counts, clearArray: false);
ulongPool.Return(Cycles, clearArray: false);
Divisors = Array.Empty<ulong>();
Exponents = Array.Empty<ulong>();
FilteredDivisors = Array.Empty<ulong>();
DivisorData = Array.Empty<GpuDivisorPartialData>();
Offsets = Array.Empty<int>();
Counts = Array.Empty<int>();
@@ -37,7 +37,7 @@ public sealed class MersenneNumberDivisorGpuTester
public bool IsDivisible(ulong exponent, in ReadOnlyGpuUInt128 divisor)
{
var gpu = GpuContextPool.RentPreferred(preferCpu: false);
var gpu = GpuContextPool.Rent();
var accelerator = gpu.Accelerator;
var kernel = GetKernel(accelerator);
var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D<byte>(1));
@@ -17,7 +17,6 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
throw new NotImplementedException($"GPU incremental scanning requires the device cycle heuristics implementation (kernel: {_kernelType}, GPU order: {_useGpuOrder}).");
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope();
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy
@@ -47,88 +46,82 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
int currentSize;
int i;
UInt128 q = UInt128.Zero;
try
while (kStart <= maxK && Volatile.Read(ref isPrime))
{
while (kStart <= maxK && Volatile.Read(ref isPrime))
remaining = maxK - kStart + UInt128.One;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize);
// Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + UInt128.One;
if (_kernelType == GpuKernelType.Pow2Mod)
{
remaining = maxK - kStart + UInt128.One;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize);
// Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + UInt128.One;
if (_kernelType == GpuKernelType.Pow2Mod)
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
pow2Kernel(
stream,
currentSize,
exponent,
twoPGpu,
(GpuUInt128)kStart,
last,
divMul,
kernelArgs,
orderBuffer.View,
smallCyclesView,
primeViews.LastOne,
primeViews.LastSeven,
primeViews.LastOnePow2,
primeViews.LastSevenPow2);
}
else
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
incKernel(
stream,
currentSize,
exponent,
twoPGpu,
(GpuUInt128)kStart,
last,
divMul,
q0m10,
q0m8,
q0m3,
q0m5,
orderBuffer.View,
smallCyclesView);
}
stream.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
if (_kernelType == GpuKernelType.Pow2Mod)
{
for (i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
}
else
{
q = twoP * kStart + UInt128.One;
for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
{
if (orders[i] != 0UL && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
}
}
}
kStart += (UInt128)currentSize;
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
pow2Kernel(
stream,
currentSize,
exponent,
twoPGpu,
(GpuUInt128)kStart,
last,
divMul,
kernelArgs,
orderBuffer.View,
smallCyclesView,
primeViews.LastOne,
primeViews.LastSeven,
primeViews.LastOnePow2,
primeViews.LastSevenPow2);
}
else
{
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
incKernel(
stream,
currentSize,
exponent,
twoPGpu,
(GpuUInt128)kStart,
last,
divMul,
q0m10,
q0m8,
q0m3,
q0m5,
orderBuffer.View,
smallCyclesView);
}
stream.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
if (_kernelType == GpuKernelType.Pow2Mod)
{
for (i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
}
else
{
q = twoP * kStart + UInt128.One;
for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
{
if (orders[i] != 0UL && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
}
}
}
kStart += (UInt128)currentSize;
}
finally
{
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
execution.Dispose();
gpuLease.Dispose();
}
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
gpuLease.Dispose();
}
}
@@ -38,16 +38,6 @@ public class MersenneNumberLucasLehmerGpuTester
// Configurable LL slice size to keep kernels short. Default 32.
public int SliceSize = 32;
public bool IsMersennePrime(ulong exponent)
{
// Default to global kernel preference for backward compatibility.
bool runOnGpu = !GpuContextPool.ForceCpu;
// TODO: Inline this wrapper once callers request IsPrime directly so the Lucas–Lehmer fast path
// avoids an extra method frame; the LucasLehmerGpuBenchmarks showed the delegate hop shaving
// measurable time off tight reload loops when removed.
return IsPrime(exponent, runOnGpu);
}
public bool IsPrime(ulong exponent, bool runOnGpu)
{
// Early rejections aligned with incremental/order sieves, but safe for small p:
@@ -78,7 +68,7 @@ public class MersenneNumberLucasLehmerGpuTester
}
else
{
var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu);
var gpu = GpuContextPool.Rent();
var accelerator = gpu.Accelerator;
var kernel = GetKernel(accelerator);
var modulus = new GpuUInt128(((UInt128)1 << (int)exponent) - 1UL); // TODO: Cache these Mersenne moduli per exponent so LL GPU runs skip rebuilding them every launch.
@@ -110,7 +100,7 @@ public class MersenneNumberLucasLehmerGpuTester
}
}
var gpu = GpuContextPool.RentPreferred(preferCpu: false);
var gpu = GpuContextPool.Rent();
var accelerator = gpu.Accelerator;
var kernel = GetBatchKernel(accelerator);
@@ -152,23 +142,17 @@ public class MersenneNumberLucasLehmerGpuTester
}
GpuUInt128[] buffer = ArrayPool<GpuUInt128>.Shared.Rent(count);
try
ComputeResidues(exponents, buffer.AsSpan(0, count));
for (int i = 0; i < count; i++)
{
ComputeResidues(exponents, buffer.AsSpan(0, count));
for (int i = 0; i < count; i++)
{
results[i] = buffer[i].IsZero;
}
}
finally
{
ArrayPool<GpuUInt128>.Shared.Return(buffer);
results[i] = buffer[i].IsZero;
}
ArrayPool<GpuUInt128>.Shared.Return(buffer);
}
private bool IsMersennePrimeNtt(ulong exponent, GpuUInt128 nttMod, GpuUInt128 primitiveRoot, bool runOnGpu)
{
var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu);
var gpu = GpuContextPool.Rent();
var accelerator = gpu.Accelerator;
var addKernel = GetAddSmallKernel(accelerator);
var subKernel = GetSubSmallKernel(accelerator);
@@ -4,19 +4,18 @@ namespace PerfectNumbers.Core.Gpu;
public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuOrder)
{
private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useGpuOrder = useGpuOrder;
private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useGpuOrder = useGpuOrder;
public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope();
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy and avoids TDR
UInt128 kStart = UInt128.One;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
byte last = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy and avoids TDR
UInt128 kStart = UInt128.One;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
byte last = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
var kernel = _kernelType switch
{
@@ -26,48 +25,41 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO
var foundBuffer = accelerator.Allocate1D<int>(1); // TODO: Replace this allocation with the pooled device buffer from GpuOrderKernelBenchmarks so repeated scans reuse the pinned staging memory instead of allocating per run.
exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
try
UInt128 remaining;
int currentSize;
int found = 0;
while (kStart <= maxK && Volatile.Read(ref isPrime))
{
UInt128 remaining;
int currentSize;
int found = 0;
while (kStart <= maxK && Volatile.Read(ref isPrime))
{
remaining = maxK - kStart + 1UL;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
// This clear is required because the result may not always be set. Remove it after the kernels are modified to always set the result.
foundBuffer.MemSetToZero();
// Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + 1UL;
// TODO: Replace these modulo recomputations with the divisor-cycle stepping deltas captured in
// MersenneDivisorCycleLengthGpuBenchmarks so each batch advances via cached residues from the
// single snapshot block instead of recalculating mod 10/8/5/3 in every iteration.
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
remaining = maxK - kStart + 1UL;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
// This clear is required because the result may not always be set. Remove it after the kernels are modified to always set the result.
foundBuffer.MemSetToZero();
// Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + 1UL;
// TODO: Replace these modulo recomputations with the divisor-cycle stepping deltas captured in
// MersenneDivisorCycleLengthGpuBenchmarks so each batch advances via cached residues from the
// single snapshot block instead of recalculating mod 10/8/5/3 in every iteration.
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
// Ensure device has small cycles table for early rejections in order kernels
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
kernel(stream, currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul,
kernelArgs, foundBuffer.View, smallCyclesView);
// Ensure device has small cycles table for early rejections in order kernels
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
kernel(stream, currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul,
kernelArgs, foundBuffer.View, smallCyclesView);
stream.Synchronize();
foundBuffer.View.CopyToCPU(ref found, 1);
if (found != 0)
{
Volatile.Write(ref isPrime, false);
break;
}
stream.Synchronize();
foundBuffer.View.CopyToCPU(ref found, 1);
if (found != 0)
{
Volatile.Write(ref isPrime, false);
break;
}
kStart += (UInt128)currentSize;
}
kStart += (UInt128)currentSize;
}
}
finally
{
foundBuffer.Dispose();
execution.Dispose();
gpuLease.Dispose();
}
}
foundBuffer.Dispose();
gpuLease.Dispose();
}
}
@@ -30,7 +30,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
var gpuLease = RentAccelerator();
var accelerator = gpuLease.Accelerator;
Monitor.Enter(gpuLease.ExecutionLock);
// Monitor.Enter(gpuLease.ExecutionLock);
var stream = accelerator.CreateStream();
var kernel = GetPow2ModKernel(accelerator);
@@ -108,7 +108,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
ReturnResources(accelerator, resources);
stream.Dispose();
Monitor.Exit(gpuLease.ExecutionLock);
// Monitor.Exit(gpuLease.ExecutionLock);
ReturnAccelerator(gpuLease);
limiter.Dispose();
}
@@ -132,7 +132,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
return lease;
}
return GpuContextPool.RentPreferred(preferCpu: !_useGpuOrder);
return GpuContextPool.Rent();
}
private void ReturnAccelerator(GpuContextPool.GpuContextLease lease) => _acceleratorPool.Add(lease);
@@ -161,7 +161,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
private void ReturnResources(Accelerator accelerator, ResidueResources resources) =>
_resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag<ResidueResources>()).Add(resources);
private sealed class ResidueResources : IDisposable
private sealed class ResidueResources
{
internal ResidueResources(Accelerator accelerator, int capacity)
{
+1 -1
View File
@@ -47,7 +47,7 @@ public static class NttGpuMath
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int>> BitReverseKernelCache = new(); // TODO: Precompute bit-reversal kernels during boot so we can remove the concurrent lookup overhead flagged in the benchmarks.
private sealed class SquareCacheEntry : IDisposable
private sealed class SquareCacheEntry
{
public MemoryBuffer1D<GpuUInt128, Stride1D.Dense> Buffer { get; }
public GpuUInt128 Root { get; }
@@ -194,21 +194,20 @@ internal static partial class PrimeOrderGpuHeuristics
try
{
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var kernel = GetPartialFactorKernel(accelerator);
SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator);
SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator);
// TODO: We should create / reallocate these buffer only once or if the new required length is bigger than capacity.
// TODO: We should create / reallocate these buffer only once or if the new required length is bigger than capacity.
var factorBuffer = accelerator.Allocate1D<ulong>(primeTargets.Length);
var exponentBuffer = accelerator.Allocate1D<int>(exponentTargets.Length);
var countBuffer = accelerator.Allocate1D<int>(1);
var remainingBuffer = accelerator.Allocate1D<ulong>(1);
var fullyFactoredBuffer = accelerator.Allocate1D<byte>(1);
// TODO: There is no need to clear these buffers because the kernel will always assign values within the required bounds.
// TODO: There is no need to clear these buffers because the kernel will always assign values within the required bounds.
// factorBuffer.MemSetToZero();
// exponentBuffer.MemSetToZero();
// countBuffer.MemSetToZero();
@@ -242,13 +241,12 @@ internal static partial class PrimeOrderGpuHeuristics
fullyFactoredBuffer.View.CopyToCPU(ref fullyFactoredFlag, 1);
fullyFactored = fullyFactoredFlag != 0;
// TODO: These buffers shouldn't be disposed but rather reused accross call with the assigned accelerator.
// TODO: These buffers shouldn't be disposed but rather reused accross call with the assigned accelerator.
factorBuffer.Dispose();
exponentBuffer.Dispose();
countBuffer.Dispose();
remainingBuffer.Dispose();
fullyFactoredBuffer.Dispose();
execution.Dispose();
lease.Dispose();
return true;
}
@@ -362,12 +360,11 @@ internal static partial class PrimeOrderGpuHeuristics
order = 0UL;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var kernel = GetOrderKernel(accelerator);
SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator);
SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator);
// TODO: These buffers should be created reallocated once and assigned to an accelerator. Callers should use their own
// thread static cache to prevent other threads from using taking accelerators with pre-allocated buffers if they don't
@@ -467,7 +464,6 @@ internal static partial class PrimeOrderGpuHeuristics
stackProductBuffer.Dispose();
resultBuffer.Dispose();
statusBuffer.Dispose();
execution.Dispose();
}
}
@@ -526,7 +522,6 @@ internal static partial class PrimeOrderGpuHeuristics
private static bool TryComputeOnGpu(ReadOnlySpan<ulong> exponents, ulong prime, in MontgomeryDivisorData divisorData, Span<ulong> results)
{
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var kernel = GetPow2ModKernel(accelerator);
@@ -551,7 +546,6 @@ internal static partial class PrimeOrderGpuHeuristics
Console.WriteLine($"Exception for {prime} and exponents: {string.Join(",", exponents.ToArray())}.");
exponentBuffer.Dispose();
remainderBuffer.Dispose();
execution.Dispose();
lease.Dispose();
throw;
}
@@ -559,7 +553,6 @@ internal static partial class PrimeOrderGpuHeuristics
remainderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), exponents.Length);
exponentBuffer.Dispose();
remainderBuffer.Dispose();
execution.Dispose();
lease.Dispose();
return true;
}
@@ -587,7 +580,6 @@ internal static partial class PrimeOrderGpuHeuristics
try
{
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var kernel = GetPow2ModWideKernel(accelerator);
@@ -628,7 +620,6 @@ internal static partial class PrimeOrderGpuHeuristics
exponentBuffer.Dispose();
remainderBuffer.Dispose();
execution.Dispose();
return true;
}
catch (Exception)
@@ -0,0 +1,717 @@
using Open.Numeric.Primes;
using System;
using ILGPU.Runtime;
using System.Runtime.CompilerServices;
using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core;
public sealed class HeuristicCombinedPrimeTester
{
static HeuristicCombinedPrimeTester()
{
HeuristicPrimeSieves.EnsureInitialized();
CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding1TwoAOneBSquares);
CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding3TwoAOneBSquares);
CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding7TwoAOneBSquares);
CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding9TwoAOneBSquares);
CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding1OneAOneBSquares);
CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding3OneAOneBSquares);
CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding7OneAOneBSquares);
CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding9OneAOneBSquares);
CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding1ThreeAOneBSquares);
CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding3ThreeAOneBSquares);
CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding7ThreeAOneBSquares);
CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding9ThreeAOneBSquares);
CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding1ThreeATwoBSquares);
CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding3ThreeATwoBSquares);
CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding7ThreeATwoBSquares);
CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding9ThreeATwoBSquares);
CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB;
CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB;
CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB;
CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB;
CombinedDivisorsEnding1Squares = CombinedDivisorsEnding1OneAOneBSquares;
CombinedDivisorsEnding3Squares = CombinedDivisorsEnding3OneAOneBSquares;
CombinedDivisorsEnding7Squares = CombinedDivisorsEnding7OneAOneBSquares;
CombinedDivisorsEnding9Squares = CombinedDivisorsEnding9OneAOneBSquares;
}
public static void EnsureInitialized()
{
// Intentionally left blank. Accessing this method forces the static constructor to run.
}
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding1Span => CombinedDivisorsEnding1;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding3Span => CombinedDivisorsEnding3;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding7Span => CombinedDivisorsEnding7;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding9Span => CombinedDivisorsEnding9;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding1SquaresSpan => CombinedDivisorsEnding1Squares;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding3SquaresSpan => CombinedDivisorsEnding3Squares;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding7SquaresSpan => CombinedDivisorsEnding7Squares;
internal static ReadOnlySpan<ulong> CombinedDivisorsEnding9SquaresSpan => CombinedDivisorsEnding9Squares;
private const ulong Wheel210 = 210UL;
private static readonly ulong[] CombinedDivisorsEnding1TwoAOneB;
private static readonly ulong[] CombinedDivisorsEnding3TwoAOneB;
private static readonly ulong[] CombinedDivisorsEnding7TwoAOneB;
private static readonly ulong[] CombinedDivisorsEnding9TwoAOneB;
private static readonly ulong[] CombinedDivisorsEnding1OneAOneB;
private static readonly ulong[] CombinedDivisorsEnding3OneAOneB;
private static readonly ulong[] CombinedDivisorsEnding7OneAOneB;
private static readonly ulong[] CombinedDivisorsEnding9OneAOneB;
private static readonly ulong[] CombinedDivisorsEnding1ThreeAOneB;
private static readonly ulong[] CombinedDivisorsEnding3ThreeAOneB;
private static readonly ulong[] CombinedDivisorsEnding7ThreeAOneB;
private static readonly ulong[] CombinedDivisorsEnding9ThreeAOneB;
private static readonly ulong[] CombinedDivisorsEnding1ThreeATwoB;
private static readonly ulong[] CombinedDivisorsEnding3ThreeATwoB;
private static readonly ulong[] CombinedDivisorsEnding7ThreeATwoB;
private static readonly ulong[] CombinedDivisorsEnding9ThreeATwoB;
private static readonly ulong[] CombinedDivisorsEnding1;
private static readonly ulong[] CombinedDivisorsEnding3;
private static readonly ulong[] CombinedDivisorsEnding7;
private static readonly ulong[] CombinedDivisorsEnding9;
private static readonly ulong[] CombinedDivisorsEnding1TwoAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding3TwoAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding7TwoAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding9TwoAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding1OneAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding3OneAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding7OneAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding9OneAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding1ThreeAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding3ThreeAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding7ThreeAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding9ThreeAOneBSquares;
private static readonly ulong[] CombinedDivisorsEnding1ThreeATwoBSquares;
private static readonly ulong[] CombinedDivisorsEnding3ThreeATwoBSquares;
private static readonly ulong[] CombinedDivisorsEnding7ThreeATwoBSquares;
private static readonly ulong[] CombinedDivisorsEnding9ThreeATwoBSquares;
private static readonly ulong[] CombinedDivisorsEnding1Squares;
private static readonly ulong[] CombinedDivisorsEnding3Squares;
private static readonly ulong[] CombinedDivisorsEnding7Squares;
private static readonly ulong[] CombinedDivisorsEnding9Squares;
private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
internal enum CombinedDivisorPattern : byte
{
TwoAOneB = 0,
OneAOneB = 1,
ThreeAOneB = 2,
ThreeATwoB = 3,
}
internal enum HeuristicDivisorGroup : byte
{
None = 0,
GroupAConstant = 1,
GroupAWheel = 2,
GroupB = 3,
}
internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue)
{
public readonly ulong Value = value;
public readonly HeuristicDivisorGroup Group = group;
public readonly byte Ending = ending;
public readonly byte PriorityIndex = priorityIndex;
public readonly ushort WheelResidue = wheelResidue;
}
internal readonly struct HeuristicDivisorPreparation(
in HeuristicDivisorCandidate candidate,
in MontgomeryDivisorData divisorData,
ulong cycleLengthHint,
bool hasCycleLengthHint)
{
public readonly HeuristicDivisorCandidate Candidate = candidate;
public readonly MontgomeryDivisorData DivisorData = divisorData;
public readonly ulong CycleLengthHint = cycleLengthHint;
public readonly bool HasCycleLengthHint = hasCycleLengthHint;
public readonly bool RequiresCycleComputation = !hasCycleLengthHint;
}
[ThreadStatic]
private static HeuristicCombinedPrimeTester? _tester;
public static HeuristicCombinedPrimeTester Exclusive => _tester ??= new();
public static bool IsPrimeCpu(ulong n)
{
byte nMod10 = (byte)n.Mod10();
ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
if (maxDivisorSquare < 9UL)
{
return EvaluateWithOpenNumericFallback(n);
}
return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10);
}
public static bool IsPrimeGpu(ulong n)
{
byte nMod10 = (byte)n.Mod10();
ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
return IsPrimeGpu(n, maxDivisorSquare, nMod10);
}
public static bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10)
{
// TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
if (maxDivisorSquare < 9UL)
{
return EvaluateWithOpenNumericFallback(n);
}
bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10);
return !compositeDetected;
}
private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10)
{
ReadOnlySpan<ulong> combinedDivisors = GetCombinedDivisors(nMod10);
ReadOnlySpan<ulong> combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
int length = combinedDivisors.Length;
for (int i = 0; i < length; i++)
{
ulong divisorSquare = combinedDivisorSquares[i];
if (divisorSquare > maxDivisorSquare)
{
break;
}
if (n % combinedDivisors[i] == 0UL)
{
return false;
}
}
return EvaluateWithOpenNumericFallback(n);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static int GetBatchSize(ulong n) => n switch
{
< 1_000 => 8,
< 10_000 => 16,
< 100_000 => 32,
< 1_000_000 => 128,
_ => 256,
};
private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
{
var limiter = GpuPrimeWorkLimiter.Acquire();
var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(nMod10, 1, PrimeTester.PrimeTesterGpuContextPool.HeuristicGpuTableScope.GroupABWithCombined);
var accelerator = gpu.Accelerator;
var kernel = gpu.HeuristicTrialDivisionKernel;
var flagView1D = gpu.HeuristicFlag.View;
var flagView = flagView1D.AsContiguous();
ulong maxDivisor = ComputeMaxDivisorLimit(maxDivisorSquare);
bool compositeDetected = false;
int divisorLength = GetCombinedDivisors(nMod10).Length;
int compositeFlag = 0;
flagView1D.CopyFromCPU(ref compositeFlag, 1);
kernel(
divisorLength,
flagView,
n,
maxDivisor,
HeuristicGpuDivisorTableKind.Combined,
gpu.HeuristicGpuTables,
gpu.HeuristicCombinedTables);
accelerator.Synchronize();
flagView1D.CopyToCPU(ref compositeFlag, 1);
compositeDetected = compositeFlag != 0;
gpu.Dispose();
limiter.Dispose();
return compositeDetected;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ComputeMaxDivisorLimit(ulong maxDivisorSquare)
{
ulong maxDivisor = (ulong)Math.Sqrt(maxDivisorSquare);
while ((maxDivisor + 1) != 0 && (maxDivisor + 1) <= maxDivisorSquare / (maxDivisor + 1))
{
maxDivisor++;
}
while (maxDivisor != 0 && maxDivisor > maxDivisorSquare / maxDivisor)
{
maxDivisor--;
}
return maxDivisor;
}
private static bool EvaluateWithOpenNumericFallback(ulong n)
{
return Prime.Numbers.IsPrime(n);
}
internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n)
{
return n;
}
private static ReadOnlySpan<ulong> GetCombinedDivisors(byte nMod10)
{
return nMod10 switch
{
1 => CombinedDivisorsEnding1,
3 => CombinedDivisorsEnding3,
7 => CombinedDivisorsEnding7,
9 => CombinedDivisorsEnding9,
_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
};
}
private static ReadOnlySpan<ulong> GetCombinedDivisorSquares(byte nMod10)
{
return nMod10 switch
{
1 => CombinedDivisorsEnding1Squares,
3 => CombinedDivisorsEnding3Squares,
7 => CombinedDivisorsEnding7Squares,
9 => CombinedDivisorsEnding9Squares,
_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
};
}
internal static ReadOnlySpan<ulong> GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern)
{
return pattern switch
{
CombinedDivisorPattern.TwoAOneB => nMod10 switch
{
1 => CombinedDivisorsEnding1TwoAOneB,
3 => CombinedDivisorsEnding3TwoAOneB,
7 => CombinedDivisorsEnding7TwoAOneB,
9 => CombinedDivisorsEnding9TwoAOneB,
_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
},
CombinedDivisorPattern.OneAOneB => nMod10 switch
{
1 => CombinedDivisorsEnding1OneAOneB,
3 => CombinedDivisorsEnding3OneAOneB,
7 => CombinedDivisorsEnding7OneAOneB,
9 => CombinedDivisorsEnding9OneAOneB,
_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
},
CombinedDivisorPattern.ThreeAOneB => nMod10 switch
{
1 => CombinedDivisorsEnding1ThreeAOneB,
3 => CombinedDivisorsEnding3ThreeAOneB,
7 => CombinedDivisorsEnding7ThreeAOneB,
9 => CombinedDivisorsEnding9ThreeAOneB,
_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
},
CombinedDivisorPattern.ThreeATwoB => nMod10 switch
{
1 => CombinedDivisorsEnding1ThreeATwoB,
3 => CombinedDivisorsEnding3ThreeATwoB,
7 => CombinedDivisorsEnding7ThreeATwoB,
9 => CombinedDivisorsEnding9ThreeATwoB,
_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
},
_ => throw new InvalidOperationException($"Unsupported combined divisor pattern: {pattern}."),
};
}
internal static ReadOnlySpan<ulong> GetCombinedDivisorSquares(byte nMod10, CombinedDivisorPattern pattern)
{
return pattern switch
{
CombinedDivisorPattern.TwoAOneB => nMod10 switch
{
1 => CombinedDivisorsEnding1TwoAOneBSquares,
3 => CombinedDivisorsEnding3TwoAOneBSquares,
7 => CombinedDivisorsEnding7TwoAOneBSquares,
9 => CombinedDivisorsEnding9TwoAOneBSquares,
_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
},
CombinedDivisorPattern.OneAOneB => nMod10 switch
{
1 => CombinedDivisorsEnding1OneAOneBSquares,
3 => CombinedDivisorsEnding3OneAOneBSquares,
7 => CombinedDivisorsEnding7OneAOneBSquares,
9 => CombinedDivisorsEnding9OneAOneBSquares,
_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
},
CombinedDivisorPattern.ThreeAOneB => nMod10 switch
{
1 => CombinedDivisorsEnding1ThreeAOneBSquares,
3 => CombinedDivisorsEnding3ThreeAOneBSquares,
7 => CombinedDivisorsEnding7ThreeAOneBSquares,
9 => CombinedDivisorsEnding9ThreeAOneBSquares,
_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
},
CombinedDivisorPattern.ThreeATwoB => nMod10 switch
{
1 => CombinedDivisorsEnding1ThreeATwoBSquares,
3 => CombinedDivisorsEnding3ThreeATwoBSquares,
7 => CombinedDivisorsEnding7ThreeATwoBSquares,
9 => CombinedDivisorsEnding9ThreeATwoBSquares,
_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
},
_ => throw new InvalidOperationException($"Unsupported combined divisor pattern for squares: {pattern}."),
};
}
private static ulong[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern, out ulong[] squares)
{
ReadOnlySpan<ulong> groupA = HeuristicPrimeSieves.GroupADivisors;
ReadOnlySpan<uint> groupB = GetGroupBDivisors(nMod10);
var combined = new List<ulong>(groupA.Length + groupB.Length);
int groupAIndex = 0;
while (groupAIndex < groupA.Length && groupA[groupAIndex] <= 19)
{
combined.Add(groupA[groupAIndex]);
groupAIndex++;
}
int groupBIndex = GetGroupBStartIndex(groupB);
int groupAInterleaveCount;
int groupBInterleaveCount;
switch (pattern)
{
case CombinedDivisorPattern.OneAOneB:
groupAInterleaveCount = 1;
groupBInterleaveCount = 1;
break;
case CombinedDivisorPattern.TwoAOneB:
groupAInterleaveCount = 2;
groupBInterleaveCount = 1;
break;
case CombinedDivisorPattern.ThreeAOneB:
groupAInterleaveCount = 3;
groupBInterleaveCount = 1;
break;
case CombinedDivisorPattern.ThreeATwoB:
groupAInterleaveCount = 3;
groupBInterleaveCount = 2;
break;
default:
groupAInterleaveCount = 2;
groupBInterleaveCount = 1;
break;
}
while (groupAIndex < groupA.Length || groupBIndex < groupB.Length)
{
int addedA = 0;
while (addedA < groupAInterleaveCount && groupAIndex < groupA.Length)
{
combined.Add(groupA[groupAIndex]);
groupAIndex++;
addedA++;
}
int addedB = 0;
while (addedB < groupBInterleaveCount && groupBIndex < groupB.Length)
{
combined.Add(groupB[groupBIndex]);
groupBIndex++;
addedB++;
}
}
ulong[] result = [.. combined];
ulong[] resultSquares = new ulong[result.Length];
for (int i = 0; i < result.Length; i++)
{
ulong divisor = result[i];
resultSquares[i] = divisor * divisor;
}
squares = resultSquares;
return result;
}
private static int GetGroupBStartIndex(ReadOnlySpan<uint> divisors)
{
int index = 0;
while (index < divisors.Length && IsGroupAPrefixValue(divisors[index]))
{
index++;
}
return index;
}
private static bool IsGroupAPrefixValue(ulong value) => value is 3UL or 7UL or 11UL or 13UL;
private static HeuristicDivisorGroup ResolveGroup(ulong value) => value switch
{
3UL or 7UL or 11UL or 13UL => HeuristicDivisorGroup.GroupAConstant,
_ when value % 10UL == 3UL => HeuristicDivisorGroup.GroupAWheel,
_ => HeuristicDivisorGroup.GroupB,
};
private static ReadOnlySpan<uint> GetGroupBDivisors(byte ending) => ending switch
{
1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree,
3 => [],
7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree,
9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree,
_ => [],
};
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
{
return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer);
}
internal struct HeuristicGroupBSequenceState
{
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor)
{
return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate)
{
ulong divisor = candidate.Value;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength);
return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static ulong ResolveHeuristicCycleLength(
ulong exponent,
in HeuristicDivisorPreparation preparation,
out bool cycleFromHint,
out bool cycleComputed,
out bool primeOrderFailed)
{
if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL)
{
cycleFromHint = true;
cycleComputed = true;
primeOrderFailed = false;
return preparation.CycleLengthHint;
}
ulong divisor = preparation.Candidate.Value;
MontgomeryDivisorData divisorData = preparation.DivisorData;
bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
divisor,
exponent,
divisorData,
out ulong computedCycle,
out bool primeOrderFailedLocal);
if (trySuccess && computedCycle != 0UL)
{
cycleFromHint = false;
cycleComputed = true;
primeOrderFailed = primeOrderFailedLocal;
return computedCycle;
}
primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL;
ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
cycleFromHint = false;
cycleComputed = resolvedCycle != 0UL;
return resolvedCycle;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength)
{
ulong[] snapshot = HeuristicSmallCycleSnapshot;
if (divisor < (ulong)snapshot.Length)
{
cycleLength = snapshot[(int)divisor];
if (cycleLength != 0UL)
{
return true;
}
}
cycleLength = 0UL;
return false;
}
internal ref struct HeuristicDivisorEnumerator
{
private readonly ulong maxDivisorSquare;
private readonly ReadOnlySpan<ulong> combinedDivisors;
private readonly ReadOnlySpan<ulong> combinedDivisorSquares;
private int index;
public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
{
this.maxDivisorSquare = maxDivisorSquare;
_ = groupBBuffer;
combinedDivisors = GetCombinedDivisors(nMod10);
combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
index = 0;
}
public bool TryGetNext(out HeuristicDivisorCandidate candidate)
{
while (index < combinedDivisors.Length)
{
int currentIndex = index;
ulong divisorSquare = combinedDivisorSquares[currentIndex];
if (divisorSquare > maxDivisorSquare)
{
index = combinedDivisors.Length;
break;
}
ulong entry = combinedDivisors[currentIndex];
index = currentIndex + 1;
HeuristicDivisorGroup group = ResolveGroup(entry);
ulong value = entry;
byte ending = (byte)(value % 10UL);
ushort residue = (ushort)(value % Wheel210);
candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue);
return true;
}
candidate = default;
return false;
}
}
internal ref struct MersenneHeuristicDivisorEnumerator
{
private readonly GpuUInt128 step;
private readonly GpuUInt128 limit;
private GpuUInt128 current;
private MersenneDivisorResidueStepper residueStepper;
private bool active;
private ulong processedCount;
private ulong lastDivisor;
public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor)
{
var stepLocal = new GpuUInt128(exponent);
stepLocal.ShiftLeft(1);
step = stepLocal;
var limitLocal = new GpuUInt128(maxDivisor);
limit = limitLocal;
var firstDivisor = stepLocal;
firstDivisor.Add(1UL);
bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0;
current = hasCandidates ? firstDivisor : GpuUInt128.Zero;
residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default;
active = hasCandidates;
processedCount = 0UL;
lastDivisor = 0UL;
}
public bool TryGetNext(out HeuristicDivisorCandidate candidate)
{
while (active)
{
ulong value = current.Low;
processedCount++;
lastDivisor = value;
bool admissible = residueStepper.IsAdmissible();
Advance();
if (admissible)
{
candidate = CreateCandidate(value);
return true;
}
}
candidate = default;
return false;
}
private void Advance()
{
GpuUInt128 next = current;
next.Add(step);
if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0)
{
active = false;
return;
}
current = next;
residueStepper.Advance();
}
public readonly ulong ProcessedCount => processedCount;
public readonly ulong LastDivisor => lastDivisor;
public readonly bool Exhausted => !active;
private static HeuristicDivisorCandidate CreateCandidate(ulong value)
{
byte ending = (byte)(value % 10UL);
HeuristicDivisorGroup group;
if (value == 3UL || value == 7UL || value == 11UL || value == 13UL)
{
group = HeuristicDivisorGroup.GroupAConstant;
}
else if (ending == 3)
{
group = HeuristicDivisorGroup.GroupAWheel;
}
else
{
group = HeuristicDivisorGroup.GroupB;
}
ushort residue = (ushort)(value % Wheel210);
return new HeuristicDivisorCandidate(value, group, ending, 0, residue);
}
}
}
+20 -19
View File
@@ -15,19 +15,19 @@ public static class HeuristicPrimeSieves
private static readonly ushort[] GroupBResiduesEnding7 = [17, 37, 47, 67, 97, 107, 127, 137, 157, 167, 187, 197];
private static readonly ushort[] GroupBResiduesEnding9 = [19, 29, 59, 79, 89, 109, 139, 149, 169, 179, 199, 209];
private static readonly int[] GroupADivisorsStorage;
private static readonly ulong[] GroupADivisorsStorage;
private static readonly ulong[] GroupADivisorSquaresStorage;
private static readonly int[] GroupBDivisorsEnding1Storage;
private static readonly ulong[] GroupBDivisorsEnding1Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding1Storage;
private static readonly int[] GroupBDivisorsEnding3Storage;
private static readonly ulong[] GroupBDivisorsEnding3Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding3Storage;
private static readonly int[] GroupBDivisorsEnding7Storage;
private static readonly ulong[] GroupBDivisorsEnding7Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding7Storage;
private static readonly int[] GroupBDivisorsEnding9Storage;
private static readonly ulong[] GroupBDivisorsEnding9Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding9Storage;
static HeuristicPrimeSieves()
@@ -46,45 +46,46 @@ public static class HeuristicPrimeSieves
// Intentionally left blank. Accessing this method forces the static constructor to run.
}
public static ReadOnlySpan<int> GroupADivisors => GroupADivisorsStorage;
public static ReadOnlySpan<ulong> GroupADivisors => GroupADivisorsStorage;
public static ReadOnlySpan<ulong> GroupADivisorSquares => GroupADivisorSquaresStorage;
public static ReadOnlySpan<int> GroupBDivisorsEnding1 => GroupBDivisorsEnding1Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding1 => GroupBDivisorsEnding1Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding1 => GroupBDivisorSquaresEnding1Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding3 => GroupBDivisorsEnding3Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding3 => GroupBDivisorsEnding3Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding3 => GroupBDivisorSquaresEnding3Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding7 => GroupBDivisorsEnding7Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding7 => GroupBDivisorsEnding7Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding7 => GroupBDivisorSquaresEnding7Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding9 => GroupBDivisorsEnding9Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding9 => GroupBDivisorsEnding9Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding9 => GroupBDivisorSquaresEnding9Storage;
private static (int[] divisors, ulong[] squares) BuildGroupADivisors()
private static (ulong[] divisors, ulong[] squares) BuildGroupADivisors()
{
var divisors = new int[SieveLength];
var divisors = new ulong[SieveLength];
var squares = new ulong[SieveLength];
int index = 0;
for (; index < GroupAConstants.Length && index < divisors.Length; index++)
{
int value = GroupAConstants[index];
ulong value = (ulong)GroupAConstants[index];
divisors[index] = value;
squares[index] = (ulong)value * (ulong)value;
squares[index] = value * value;
}
int candidate = 23;
int incrementIndex = 0;
while (index < divisors.Length)
{
divisors[index] = candidate;
squares[index] = (ulong)candidate * (ulong)candidate;
ulong value = (ulong)candidate;
divisors[index] = value;
squares[index] = value * value;
candidate += GroupAIncrementPattern[incrementIndex];
incrementIndex ^= 1;
@@ -94,9 +95,9 @@ public static class HeuristicPrimeSieves
return (divisors, squares);
}
private static (int[] divisors, ulong[] squares) BuildGroupBResidueSieve(ReadOnlySpan<ushort> residues)
private static (ulong[] divisors, ulong[] squares) BuildGroupBResidueSieve(ReadOnlySpan<ushort> residues)
{
var divisors = new int[SieveLength];
var divisors = new ulong[SieveLength];
var squares = new ulong[SieveLength];
int index = 0;
@@ -112,7 +113,7 @@ public static class HeuristicPrimeSieves
continue;
}
divisors[index] = (int)candidate;
divisors[index] = candidate;
squares[index] = candidate * candidate;
index++;
}
File diff suppressed because it is too large. Load diff
@@ -16,8 +16,10 @@ public interface IMersenneNumberDivisorByDivisorTester
void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues);
public interface IDivisorScanSession : IDisposable
public interface IDivisorScanSession
{
void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan<ulong> primes, Span<byte> hits);
void Dispose();
}
}
+82 -76
View File
@@ -4,6 +4,7 @@ using System.Numerics;
using System.Runtime.CompilerServices;
using System.Security.Cryptography;
using System.Text;
using System.Runtime.InteropServices;
using PerfectNumbers.Core.Gpu;
using ILGPU;
using ILGPU.Runtime;
@@ -136,20 +137,31 @@ public class MersenneDivisorCycles
{
// Binary pairs: (ulong divisor, ulong cycle)
using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true);
Stream stream = reader.BaseStream;
byte[] buffer = new byte[16];
while (true)
{
ulong d, c;
try
int read = 0;
Span<byte> span = buffer;
while (read < span.Length)
{
d = reader.ReadUInt64();
c = reader.ReadUInt64();
}
catch (EndOfStreamException)
{
yield break;
int chunk = stream.Read(span[read..]);
if (chunk == 0)
{
if (read == 0)
{
yield break;
}
throw new EndOfStreamException("Incomplete divisor-cycle entry.");
}
read += chunk;
}
yield return (d, c);
ulong divisor = BinaryPrimitives.ReadUInt64LittleEndian(span[..8]);
ulong cycle = BinaryPrimitives.ReadUInt64LittleEndian(span[8..]);
yield return (divisor, cycle);
}
}
@@ -379,7 +391,8 @@ public class MersenneDivisorCycles
for (int i = 0; i < smallLength; i++)
{
ulong prime = smallPrimes[i];
if (smallPrimesSquared[i] > remaining)
ulong primeSquared = smallPrimesSquared[i];
if (primeSquared > remaining)
{
break;
}
@@ -398,7 +411,9 @@ public class MersenneDivisorCycles
// Atomic.Add(ref _tryFactorIntoCountsInternalHits, 1UL);
// Console.WriteLine($"MersenneDivisorCycles.TryFactorIntoCountsInternal hits {Volatile.Read(ref _tryFactorIntoCountsInternalHits)}");
if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None))
// HeuristicPrimeTester primeTester = _primeTester ??= new();
if (PrimeTester.IsPrime(remaining))
// if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None))
{
AddFactor(counts, remaining);
return true;
@@ -438,26 +453,14 @@ public class MersenneDivisorCycles
private static void AddFactor(Dictionary<ulong, int> counts, ulong factor)
{
if (counts.TryGetValue(factor, out int existing))
{
counts[factor] = existing + 1;
}
else
{
counts[factor] = 1;
}
ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists);
entry = exists ? entry + 1 : 1;
}
private static void AddFactor(Dictionary<ulong, int> counts, ulong factor, int multiplicity)
{
if (counts.TryGetValue(factor, out int existing))
{
counts[factor] = existing + multiplicity;
}
else
{
counts[factor] = multiplicity;
}
ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists);
entry = exists ? entry + multiplicity : multiplicity;
}
private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
@@ -535,8 +538,8 @@ public class MersenneDivisorCycles
private static void ProcessReduceOrderPrime(Span<ulong> candidateBuffer, Span<bool> evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper)
{
ulong working = order;
int actual = 0;
while (actual < multiplicity)
int writeIndex = multiplicity;
while (writeIndex > 0)
{
if (working < prime)
{
@@ -549,26 +552,20 @@ public class MersenneDivisorCycles
break;
}
candidateBuffer[actual] = reduced;
writeIndex--;
candidateBuffer[writeIndex] = reduced;
working = reduced;
actual++;
}
int actual = multiplicity - writeIndex;
if (actual == 0)
{
return;
}
Span<ulong> candidates = candidateBuffer.Slice(0, actual);
Span<ulong> candidates = candidateBuffer.Slice(writeIndex, actual);
Span<bool> evaluations = evaluationBuffer.Slice(0, actual);
for (int left = 0, right = actual - 1; left < right; left++, right--)
{
ulong tmp = candidates[left];
candidates[left] = candidates[right];
candidates[right] = tmp;
}
stepper.Reset();
evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]);
@@ -598,24 +595,28 @@ public class MersenneDivisorCycles
// return 2UL;
// }
ulong modulus = n;
ulong nMinus1 = modulus - 1UL;
ulong nMinus2 = modulus - 2UL;
while (true)
{
ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL;
ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL;
ulong c = (NextRandomUInt64() % nMinus1) + 1UL;
ulong x = (NextRandomUInt64() % nMinus2) + 2UL;
ulong y = x;
ulong d = 1UL;
while (d == 1UL)
{
x = AdvancePolynomial(x, c, n);
y = AdvancePolynomial(y, c, n);
y = AdvancePolynomial(y, c, n);
x = AdvancePolynomial(x, c, modulus);
y = AdvancePolynomial(y, c, modulus);
y = AdvancePolynomial(y, c, modulus);
ulong diff = x > y ? x - y : y - x;
d = BinaryGcd(diff, n);
d = BinaryGcd(diff, modulus);
}
if (d != n)
if (d != modulus)
{
return d;
}
@@ -916,19 +917,22 @@ public class MersenneDivisorCycles
return CalculateCycleLengthFallback(divisor);
}
internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false)
{
if ((divisor & (divisor - 1UL)) == 0UL)
{
cycleLength = 1UL;
return true;
}
// [ThreadStatic]
// private static HeuristicPrimeTester? _primeTester;
if (divisor <= 3UL)
{
cycleLength = CalculateCycleLengthFallback(divisor);
return true;
}
internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false)
{
if ((divisor & (divisor - 1UL)) == 0UL)
{
cycleLength = 1UL;
return true;
}
if (divisor <= 3UL)
{
cycleLength = CalculateCycleLengthFallback(divisor);
return true;
}
// if (!skipPrimeOrderHeuristic)
// {
@@ -936,25 +940,27 @@ public class MersenneDivisorCycles
// Console.WriteLine($"MersenneDivisorCycles.CalculateCycleLengthFallback hits {Volatile.Read(ref _tryCalculateCycleLengthHeuristicHits)}");
// }
// if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None))
if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeCpu(divisor, CancellationToken.None))
{
ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (computedOrder != 0UL)
{
cycleLength = computedOrder;
return true;
}
}
// HeuristicPrimeTester primeTester = _primeTester ??= new();
// if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None))
// if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor))
if (!skipPrimeOrderHeuristic && PrimeTester.IsPrime(divisor))
{
ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (computedOrder != 0UL)
{
cycleLength = computedOrder;
return true;
}
}
cycleLength = 0UL;
return false;
}
cycleLength = 0UL;
return false;
}
private static ulong CalculateCycleLengthFallback(ulong divisor)
{
@@ -1,5 +1,6 @@
using System.Buffers;
using System.Collections.Concurrent;
using System.Threading;
using System.Threading.Tasks;
using System.Globalization;
using System.Runtime.InteropServices;
@@ -93,9 +94,15 @@ public static class MersenneNumberDivisorByDivisorTester
}
else
{
while (startIndex > 0 && candidates[startIndex - 1] >= startPrime)
while (startIndex > 0)
{
startIndex--;
int previousIndex = startIndex - 1;
if (candidates[previousIndex] < startPrime)
{
break;
}
startIndex = previousIndex;
}
}
}
@@ -161,14 +168,10 @@ public static class MersenneNumberDivisorByDivisorTester
if (primesToTest.Count == 0)
{
if (applyStartPrime)
{
Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor.");
}
else
{
Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering.");
}
string message = applyStartPrime
? $"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."
: "No prime candidates remain for --mersenne=bydivisor after filtering.";
Console.WriteLine(message);
return;
}
@@ -239,19 +242,20 @@ public static class MersenneNumberDivisorByDivisorTester
void ProcessPrime(ulong prime)
{
Console.WriteLine($"Task started {prime}");
Console.WriteLine($"Processing {prime}");
bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted);
if (!isPrime)
{
markComposite();
printResult(prime, true, true, false);
Console.WriteLine($"Finished processing {prime}");
return;
}
clearComposite();
printResult(prime, true, divisorsExhausted, true);
Console.WriteLine($"Task finished {prime}");
Console.WriteLine($"Finished processing {prime}");
}
if (workerCount == 1)
@@ -263,12 +267,7 @@ public static class MersenneNumberDivisorByDivisorTester
}
else
{
var scheduler = UnboundedTaskScheduler.Instance;
ParallelOptions options = new()
{
MaxDegreeOfParallelism = workerCount,
TaskScheduler = scheduler
};
TaskScheduler scheduler = UnboundedTaskScheduler.Instance;
int totalCount = filteredPrimes.Count;
int partitionSize = chunkSize < 1 ? 1 : chunkSize;
@@ -277,13 +276,41 @@ public static class MersenneNumberDivisorByDivisorTester
partitionSize = totalCount;
}
Parallel.ForEach(Partitioner.Create(0, totalCount, partitionSize), options, range =>
int taskCount = (totalCount + partitionSize - 1) / partitionSize;
Task[] tasks = new Task[taskCount];
var startGate = new ManualResetEventSlim(initialState: false);
int taskIndex = 0;
for (int start = 0; start < totalCount; start += partitionSize)
{
for (int index = range.Item1; index < range.Item2; index++)
int rangeStart = start;
int rangeEnd = rangeStart + partitionSize;
if (rangeEnd > totalCount)
{
ProcessPrime(filteredPrimes[index]);
rangeEnd = totalCount;
}
});
tasks[taskIndex++] = Task.Factory.StartNew(
() =>
{
startGate.Wait();
Console.WriteLine($"Task started for range {rangeStart}");
for (int index = rangeStart; index < rangeEnd; index++)
{
ProcessPrime(filteredPrimes[index]);
}
Console.WriteLine($"Task finished for range {rangeStart}");
},
CancellationToken.None,
TaskCreationOptions.DenyChildAttach,
scheduler);
}
startGate.Set();
Task.WaitAll(tasks);
startGate.Dispose();
}
}
}
+40 -47
View File
@@ -147,61 +147,54 @@ public sealed class MersenneNumberTester(
UInt128[] qs = qsBuffer;
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope();
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
var orderKernel = gpuLease.OrderKernel!;
try
// Guard long-running kernels by chunking the warm-up across batches.
// Reuse the same ScanBatchSize knob used for scanning.
int batchSize = GpuConstants.ScanBatchSize;
ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL;
int offset = 0;
while (offset < idx)
{
// Guard long-running kernels by chunking the warm-up across batches.
// Reuse the same ScanBatchSize knob used for scanning.
int batchSize = GpuConstants.ScanBatchSize;
ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL;
int offset = 0;
while (offset < idx)
int count = Math.Min(batchSize, idx - offset);
var qBuffer = accelerator.Allocate1D<PerfectNumbers.Core.Gpu.GpuUInt128>(count);
var orderBuffer = accelerator.Allocate1D<ulong>(count);
// Convert to device-friendly GpuUInt128 on the fly
var tmp = System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Rent(count);
for (int i = 0; i < count; i++)
{
int count = Math.Min(batchSize, idx - offset);
var qBuffer = accelerator.Allocate1D<PerfectNumbers.Core.Gpu.GpuUInt128>(count);
var orderBuffer = accelerator.Allocate1D<ulong>(count);
// Convert to device-friendly GpuUInt128 on the fly
var tmp = System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Rent(count);
for (int i = 0; i < count; i++)
{
tmp[i] = (PerfectNumbers.Core.Gpu.GpuUInt128)qs[offset + i];
}
// Copy only the populated elements to the device buffer to avoid overrun on pooled arrays.
qBuffer.View.CopyFromCPU(ref tmp[0], count);
System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Return(tmp);
orderKernel(stream, count, exponent, divMul, qBuffer.View, orderBuffer.View);
stream.Synchronize();
ulong[] orders = orderBuffer.GetAsArray1D();
for (int i = 0; i < count; i++)
{
ulong order = orders[i];
if (order == 0UL)
{
order = qs[offset + i].CalculateOrder();
}
if (useOrderCache)
{
OrderCache[qs[offset + i]] = order;
}
}
offset += count;
orderBuffer.Dispose();
qBuffer.Dispose();
tmp[i] = (PerfectNumbers.Core.Gpu.GpuUInt128)qs[offset + i];
}
// Copy only the populated elements to the device buffer to avoid overrun on pooled arrays.
qBuffer.View.CopyFromCPU(ref tmp[0], count);
System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Return(tmp);
orderKernel(stream, count, exponent, divMul, qBuffer.View, orderBuffer.View);
stream.Synchronize();
ulong[] orders = orderBuffer.GetAsArray1D();
for (int i = 0; i < count; i++)
{
ulong order = orders[i];
if (order == 0UL)
{
order = qs[offset + i].CalculateOrder();
}
if (useOrderCache)
{
OrderCache[qs[offset + i]] = order;
}
}
offset += count;
orderBuffer.Dispose();
qBuffer.Dispose();
}
finally
{
execution.Dispose();
gpuLease.Dispose();
ArrayPool<UInt128>.Shared.Return(qs);
}
gpuLease.Dispose();
ArrayPool<UInt128>.Shared.Return(qs);
}
public bool IsMersennePrime(ulong exponent)
@@ -14,7 +14,7 @@ public static class MersennePrimeFactorTester
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeFactor(ulong p, ulong q, CancellationToken ct)
{
if (q < 2UL || (q & 1UL) == 0)
if (q < 2UL)
{
return false;
}
@@ -34,7 +34,7 @@ public static class MersennePrimeFactorTester
// Atomic.Add(ref _isPrimeFactorHits, 1UL);
// Console.WriteLine($"MersennePrimeFactorTester.IsPrimeFactor hits {Volatile.Read(ref _isPrimeFactorHits)}");
if (!PrimeTester.IsPrimeCpu(q, ct))
if (!PrimeTester.IsPrime(q))
{
return false;
}
@@ -247,7 +247,7 @@ public static class MersennePrimeFactorTester
// Console.WriteLine($"MersennePrimeFactorTester.Factor64 hits {Volatile.Read(ref _factor64Hits)}");
// if (HeuristicPrimeTester.Exclusive.IsPrime(m, ct))
if (PrimeTester.IsPrimeCpu(m, ct))
if (PrimeTester.IsPrime(m))
{
dict[m] = 1;
continue;
+1 -1
View File
@@ -34,7 +34,7 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
+150 -222
View File
@@ -299,22 +299,16 @@ internal static partial class PrimeOrderCalculator
Span<ulong> candidates = candidateBuffer[..actual];
Span<bool> evaluations = evaluationBuffer[..actual];
for (int left = 0, right = actual - 1; left < right; left++, right--)
{
ulong tmp = candidates[left];
candidates[left] = candidates[right];
candidates[right] = tmp;
}
stepper.Reset();
evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]);
for (int j = 1; j < actual; j++)
int last = actual - 1;
evaluations[last] = stepper.InitializeCpuIsUnity(candidates[last]);
for (int j = last - 1; j >= 0; j--)
{
evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]);
}
for (int j = actual - 1; j >= 0; j--)
for (int j = 0; j < actual; j++)
{
if (!evaluations[j])
{
@@ -464,16 +458,16 @@ internal static partial class PrimeOrderCalculator
}
Span<ulong> candidates = buffer[..actual];
candidates.Reverse();
stepper.Reset();
if (stepper.InitializeCpuIsUnity(candidates[0]))
int last = actual - 1;
if (stepper.InitializeCpuIsUnity(candidates[last]))
{
return true;
}
for (int j = 1; j < actual; j++)
for (int j = last - 1; j >= 0; j--)
{
if (stepper.ComputeNextIsUnity(candidates[j]))
{
@@ -694,13 +688,16 @@ internal static partial class PrimeOrderCalculator
private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates)
{
ulong previous = previousOrder ?? 0UL;
int previousGroup = previousOrder.HasValue ? GetGroup(previousOrder.Value, prime) : 1;
bool hasPrevious = previousOrder.HasValue;
ulong threshold1 = prime >> 3;
ulong threshold2 = prime >> 2;
ulong threshold3 = (prime * 3UL) >> 3;
int previousGroup = hasPrevious ? GetGroup(previous, threshold1, threshold2, threshold3) : 1;
candidates.Sort((x, y) =>
{
CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious);
CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious);
CandidateKey keyX = BuildKey(x, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3);
CandidateKey keyY = BuildKey(y, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3);
int compare = keyX.Primary.CompareTo(keyY.Primary);
if (compare != 0)
{
@@ -717,9 +714,9 @@ internal static partial class PrimeOrderCalculator
});
}
private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious)
private static CandidateKey BuildKey(ulong value, ulong previous, int previousGroup, bool hasPrevious, ulong threshold1, ulong threshold2, ulong threshold3)
{
int group = GetGroup(value, prime);
int group = GetGroup(value, threshold1, threshold2, threshold3);
if (group == 0)
{
return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue);
@@ -769,21 +766,18 @@ internal static partial class PrimeOrderCalculator
return groupOffset + (isGe ? 0 : 1);
}
private static int GetGroup(ulong value, ulong prime)
private static int GetGroup(ulong value, ulong threshold1, ulong threshold2, ulong threshold3)
{
ulong threshold1 = prime >> 3;
if (value <= threshold1)
{
return 1;
}
ulong threshold2 = prime >> 2;
if (value <= threshold2)
{
return 2;
}
ulong threshold3 = (prime * 3UL) >> 3;
if (value <= threshold3)
{
return 3;
@@ -828,9 +822,11 @@ internal static partial class PrimeOrderCalculator
}
FactorEntry factor = factors[index];
int factorExponent = factor.Exponent;
ulong primeFactor = factor.Value;
ulong contribution = 1UL;
for (int exponent = 0; exponent <= factor.Exponent; exponent++)
ulong contributionLimit = factorExponent == 0 ? 0UL : order / primeFactor;
for (int exponent = 0; exponent <= factorExponent; exponent++)
{
ulong nextDivisor = divisorProduct * contribution;
if (nextDivisor > order)
@@ -844,12 +840,12 @@ internal static partial class PrimeOrderCalculator
return;
}
if (exponent == factor.Exponent)
if (exponent == factorExponent)
{
break;
}
if (contribution > order / primeFactor)
if (contribution > contributionLimit)
{
break;
}
@@ -973,24 +969,45 @@ internal static partial class PrimeOrderCalculator
}
Span<ulong> candidates = buffer[..actual];
candidates.Reverse();
stepper.Reset();
int last = actual - 1;
if (powUsed >= powBudget && powBudget > 0)
bool enforceBudget = powBudget > 0;
if (!enforceBudget)
{
powUsed++;
if (stepper.InitializeCpuIsUnity(candidates[last]))
{
return true;
}
for (int j = last - 1; j >= 0; j--)
{
powUsed++;
if (stepper.ComputeNextIsUnity(candidates[j]))
{
return true;
}
}
return false;
}
if (powUsed >= powBudget)
{
return true;
}
powUsed++;
if (stepper.InitializeCpuIsUnity(candidates[0]))
if (stepper.InitializeCpuIsUnity(candidates[last]))
{
return true;
}
for (int j = 1; j < actual; j++)
for (int j = last - 1; j >= 0; j--)
{
if (powUsed >= powBudget && powBudget > 0)
if (powUsed >= powBudget)
{
return true;
}
@@ -1040,13 +1057,13 @@ internal static partial class PrimeOrderCalculator
// Span<int> exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount);
Span<ulong> primeSlots = stackalloc ulong[FactorSlotCount];
Span<int> exponentSlots = stackalloc int[FactorSlotCount];
ulong[]? primeSlotArray = null;
int[]? exponentSlotArray = null;
// We don't need to worry about leftovers, because we always use indexes within the calculated counts
// primeSlots.Clear();
// exponentSlots.Clear();
int factorCount = 0;
Dictionary<ulong, int>? counts = null;
bool useDictionary = false;
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
ulong remaining = value;
@@ -1061,45 +1078,20 @@ internal static partial class PrimeOrderCalculator
Stack<ulong>? compositeStack = null;
PartialFactorResult result;
PrimeTester primeTester = PrimeTester.Exclusive;
CancellationToken ct = CancellationToken.None;
// HeuristicPrimeTester? tester = _tester ??= new();
if (!gpuFactored)
{
// We don't need to worry about leftovers, because we always use indexes within the calculated counts
// primeSlots.Clear();
// exponentSlots.Clear();
counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(FactorSlotCount, 8));
// counts.Clear();
// bool gpuPopulated = TryPopulateSmallPrimeFactorsGpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining);
// if (!gpuPopulated)
// {
// throw new InvalidOperationException($"GPU didn't populate factors for {value}");
// counts.Clear();
remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
// }
int dictionaryCount = counts.Count;
if (dictionaryCount <= FactorSlotCount)
{
int copyIndex = 0;
foreach (KeyValuePair<ulong, int> entry in counts)
{
primeSlots[copyIndex] = entry.Key;
exponentSlots[copyIndex] = entry.Value;
copyIndex++;
}
factorCount = copyIndex;
counts.Clear();
ThreadStaticPools.ReturnUlongIntDictionary(counts);
counts = null;
}
else
{
useDictionary = true;
factorCount = dictionaryCount;
}
PopulateSmallPrimeFactorsCpu(
value,
limit,
ref primeSlots,
ref exponentSlots,
ref primeSlotArray,
ref exponentSlotArray,
out factorCount,
out remaining);
}
@@ -1133,12 +1125,13 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorHits, 1UL);
// Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}");
bool isPrime = PrimeTester.IsPrimeGpu(composite);
// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
bool isPrime = PrimeTester.IsPrime(composite);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
if (isPrime)
{
AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1);
continue;
}
@@ -1200,7 +1193,8 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorPendingHits, 1UL);
// Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}");
bool isPrime = PrimeTester.IsPrimeGpu(composite);
// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
bool isPrime = PrimeTester.IsPrime(composite);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
entry = entry.WithPrimality(isPrime);
@@ -1209,7 +1203,7 @@ internal static partial class PrimeOrderCalculator
if (entry.IsPrime)
{
AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1);
}
// composite is never smaller on the execution path
// else if (composite > 1UL)
@@ -1237,39 +1231,22 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
cofactorIsPrime = PrimeTester.IsPrimeCpu(cofactor, CancellationToken.None);
cofactorIsPrime = PrimeTester.IsPrime(cofactor);
// cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
}
ArrayPool<FactorEntry> pool = ThreadStaticPools.FactorEntryPool;
bool temp;
if (useDictionary)
if (factorCount == 0)
{
temp = counts is null || counts.Count == 0;
if (temp && cofactor == value)
if (cofactor == value)
{
result = PartialFactorResult.Rent(null, cofactor, false, 0, cofactorIsPrime);
goto ReturnResult;
}
if (temp)
{
result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0, cofactorIsPrime);
goto ReturnResult;
}
FactorEntry[] factors = pool.Rent(counts!.Count);
index = 0;
foreach (KeyValuePair<ulong, int> entry in counts)
{
factors[index] = new FactorEntry(entry.Key, entry.Value, true);
index++;
}
Array.Sort(factors, static (a, b) => a.Value.CompareTo(b.Value));
temp = cofactor == 1UL;
result = PartialFactorResult.Rent(factors, cofactor, temp, index, cofactorIsPrime);
result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0, cofactorIsPrime);
goto ReturnResult;
}
@@ -1309,14 +1286,10 @@ internal static partial class PrimeOrderCalculator
result = PartialFactorResult.Rent(array, cofactor, temp, arrayIndex, cofactorIsPrime);
ReturnResult:
// ThreadStaticPools.UlongPool.Return(primeSlotsArray);
// ThreadStaticPools.IntPool.Return(exponentSlotsArray);
if (counts is not null)
if (primeSlotArray is not null)
{
counts.Clear();
ThreadStaticPools.ReturnUlongIntDictionary(counts);
ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false);
ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false);
}
if (compositeStack is not null)
@@ -1330,35 +1303,28 @@ internal static partial class PrimeOrderCalculator
return result;
}
private static bool TryPopulateSmallPrimeFactorsCpu(
ulong value,
uint limit,
Span<ulong> primeTargets,
Span<int> exponentTargets,
out int factorCount,
out ulong remaining)
private static void PopulateSmallPrimeFactorsCpu(
ulong value,
uint limit,
ref Span<ulong> primeTargets,
ref Span<int> exponentTargets,
ref ulong[]? primeArray,
ref int[]? exponentArray,
out int factorCount,
out ulong remaining)
{
factorCount = 0;
remaining = value;
ulong remainingLocal = value;
// value will never <= 1 in production code
// if (value <= 1UL)
// {
// return true;
// }
int capacity = Math.Min(primeTargets.Length, exponentTargets.Length);
// capacity will never equal 0 in production code
// if (capacity == 0)
// {
// return false;
// }
if (value <= 1UL)
{
remaining = remainingLocal;
return;
}
uint[] primes = PrimesGenerator.SmallPrimes;
ulong[] squares = PrimesGenerator.SmallPrimesPow2;
int primeCount = primes.Length;
ulong remainingLocal = value;
uint effectiveLimit = limit == 0 ? uint.MaxValue : limit;
for (int i = 0; i < primeCount && remainingLocal > 1UL; i++)
@@ -1377,13 +1343,13 @@ internal static partial class PrimeOrderCalculator
break;
}
ulong primeValue = primeCandidate;
// primeCandidate will never equal 0 in production code
// if (primeCandidate == 0UL)
// {
// continue;
// }
ulong primeValue = primeCandidate;
if ((remainingLocal % primeValue) != 0UL)
{
continue;
@@ -1395,10 +1361,13 @@ internal static partial class PrimeOrderCalculator
continue;
}
if (factorCount >= capacity)
{
throw new InvalidOperationException($"Capacity is smaller than factor count");
}
EnsureCollectorCapacity(
ref primeTargets,
ref exponentTargets,
ref primeArray,
ref exponentArray,
factorCount,
factorCount + 1);
primeTargets[factorCount] = primeValue;
exponentTargets[factorCount] = exponent;
@@ -1406,7 +1375,6 @@ internal static partial class PrimeOrderCalculator
}
remaining = remainingLocal;
return true;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -1474,46 +1442,6 @@ internal static partial class PrimeOrderCalculator
return exponent;
}
private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary<ulong, int> counts)
{
ulong remaining = value;
uint[] primes = PrimesGenerator.SmallPrimes;
ulong[] squares = PrimesGenerator.SmallPrimesPow2;
int primeCount = primes.Length;
for (int i = 0; i < primeCount; i++)
{
uint primeCandidate = primes[i];
if (primeCandidate > limit)
{
break;
}
if (squares[i] > remaining)
{
break;
}
if ((remaining % primeCandidate) != 0UL)
{
continue;
}
int exponent = 0;
do
{
// TODO: Implement DivRem like solution or residue stepper to re-use previous iteration results in new calculation. Add / use method, likely to ULongExtensions
remaining /= primeCandidate;
exponent++;
}
while ((remaining % primeCandidate) == 0UL);
counts[primeCandidate] = exponent;
}
return remaining;
}
private static long CreateDeadlineTimestamp(int milliseconds)
{
if (milliseconds <= 0)
@@ -1636,11 +1564,18 @@ internal static partial class PrimeOrderCalculator
}
}
private static bool TryAppendFactor(Span<ulong> primes, Span<int> exponents, ref int count, ulong prime, int exponent)
private static void AddFactorToCollector(
ref Span<ulong> primes,
ref Span<int> exponents,
ref ulong[]? primeArray,
ref int[]? exponentArray,
ref int count,
ulong prime,
int exponent)
{
if (exponent <= 0)
{
return true;
return;
}
for (int i = 0; i < count; i++)
@@ -1648,70 +1583,60 @@ internal static partial class PrimeOrderCalculator
if (primes[i] == prime)
{
exponents[i] += exponent;
return true;
return;
}
}
if (count >= primes.Length)
{
return false;
}
EnsureCollectorCapacity(
ref primes,
ref exponents,
ref primeArray,
ref exponentArray,
count,
count + 1);
primes[count] = prime;
exponents[count] = exponent;
count++;
return true;
}
private static void AddFactorToCollector(
ref bool useDictionary,
ref Dictionary<ulong, int>? counts,
Span<ulong> primes,
Span<int> exponents,
ref int count,
ulong prime,
int exponent)
private static void EnsureCollectorCapacity(
ref Span<ulong> primes,
ref Span<int> exponents,
ref ulong[]? primeArray,
ref int[]? exponentArray,
int count,
int requiredCapacity)
{
if (exponent <= 0)
if (requiredCapacity <= primes.Length)
{
return;
}
if (useDictionary)
int newCapacity = primes.Length << 1;
if (newCapacity < requiredCapacity)
{
AddFactor(counts!, prime, exponent);
return;
newCapacity = requiredCapacity;
}
if (TryAppendFactor(primes, exponents, ref count, prime, exponent))
ulong[] newPrimeArray = ThreadStaticPools.UlongPool.Rent(newCapacity);
int[] newExponentArray = ThreadStaticPools.IntPool.Rent(newCapacity);
Span<ulong> newPrimeSpan = newPrimeArray.AsSpan();
Span<int> newExponentSpan = newExponentArray.AsSpan();
primes.Slice(0, count).CopyTo(newPrimeSpan);
exponents.Slice(0, count).CopyTo(newExponentSpan);
if (primeArray is not null)
{
return;
ThreadStaticPools.UlongPool.Return(primeArray, clearArray: false);
ThreadStaticPools.IntPool.Return(exponentArray!, clearArray: false);
}
counts ??= ThreadStaticPools.RentUlongIntDictionary(Math.Max(count, 8));
CopyFactorsToDictionary(primes, exponents, count, counts);
count = 0;
useDictionary = true;
AddFactor(counts, prime, exponent);
}
private static void CopyFactorsToDictionary(
in ReadOnlySpan<ulong> primes,
in ReadOnlySpan<int> exponents,
int count,
Dictionary<ulong, int> target)
{
for (int i = 0; i < count; i++)
{
ulong prime = primes[i];
int exponent = exponents[i];
if (prime == 0UL || exponent == 0)
{
continue;
}
target[prime] = exponent;
}
primeArray = newPrimeArray;
exponentArray = newExponentArray;
primes = newPrimeSpan;
exponents = newExponentSpan;
}
private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData)
@@ -1758,9 +1683,11 @@ internal static partial class PrimeOrderCalculator
{
FactorCompletely(value, counts, knownComposite: false);
}
// private static ulong _factorCompletelyHits;
// private static ulong _factorCompletelyPollardRhoHits;
// [ThreadStatic]
// private static HeuristicPrimeTester? _tester;
private static void FactorCompletely(ulong value, Dictionary<ulong, int> counts, bool knownComposite)
{
@@ -1771,11 +1698,11 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _factorCompletelyHits, 1UL);
// Console.WriteLine($"Factor completely hits {Volatile.Read(ref _factorCompletelyHits)}");
// PrimeTester primeTester = PrimeTester.Exclusive;
// HeuristicPrimeTester tester = _tester ??= new();
// bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
// bool isPrime = PrimeTester.IsPrimeCpu(value, CancellationToken.None);
bool isPrime = PrimeTester.IsPrimeCpu(value, CancellationToken.None);
bool isPrime = PrimeTester.IsPrime(value);
// bool isPrime = PrimeTester.IsPrimeGpu(value);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
if (!knownComposite && isPrime)
{
@@ -1790,8 +1717,8 @@ internal static partial class PrimeOrderCalculator
// Console.WriteLine($"Factor completely after PollardRho hits {Volatile.Read(ref _factorCompletelyPollardRhoHits)}");
// isPrime = primeTester.HeuristicIsPrimeGpu(factor);
// isPrime = PrimeTester.IsPrimeCpu(factor, CancellationToken.None);
isPrime = PrimeTester.IsPrimeCpu(factor, CancellationToken.None);
isPrime = PrimeTester.IsPrime(factor);
// isPrime = PrimeTester.IsPrimeGpu(factor);
// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
if (isPrime)
{
@@ -1827,7 +1754,8 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
isPrime = PrimeTester.IsPrimeCpu(quotient, CancellationToken.None);
isPrime = PrimeTester.IsPrime(quotient);
// isPrime = PrimeTester.IsPrimeGpu(quotient);
// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
}
+48 -71
View File
@@ -9,86 +9,63 @@ namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
private const int GpuSmallPrimeFactorSlots = 64;
private const int GpuSmallPrimeFactorSlots = 64;
private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out int factorCount, out ulong remaining)
{
var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots);
var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots);
Span<ulong> primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
Span<int> exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
// primeBuffer.Clear();
// exponentBuffer.Clear();
remaining = value;
private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out int factorCount, out ulong remaining)
{
var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots);
var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots);
Span<ulong> primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
Span<int> exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
remaining = value;
try
{
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
try
{
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator);
SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots);
// scratch.Clear();
SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator);
SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots);
var kernel = lease.SmallPrimeFactorKernel;
kernel(
stream,
1,
value,
limit,
tables.PrimesView,
tables.SquaresView,
tables.Count,
scratch.PrimeSlots.View,
scratch.ExponentSlots.View,
scratch.CountSlot.View,
scratch.RemainingSlot.View);
var kernel = lease.SmallPrimeFactorKernel;
kernel(
stream,
1,
value,
limit,
tables.PrimesView,
tables.SquaresView,
tables.Count,
scratch.PrimeSlots.View,
scratch.ExponentSlots.View,
scratch.CountSlot.View,
scratch.RemainingSlot.View);
stream.Synchronize();
stream.Synchronize();
factorCount = 0;
scratch.CountSlot.View.CopyToCPU(ref factorCount, 1);
factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots);
factorCount = 0;
scratch.CountSlot.View.CopyToCPU(ref factorCount, 1);
factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots);
if (factorCount > 0)
{
scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount);
scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount);
}
if (factorCount > 0)
{
scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount);
scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount);
}
scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1);
scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1);
for (int i = 0; i < factorCount; i++)
{
ulong primeValue = primeBuffer[i];
int exponent = exponentBuffer[i];
// This will never happen in production code
// if (primeValue == 0UL || exponent == 0)
// {
// continue;
// }
for (int i = 0; i < factorCount; i++)
{
ulong primeValue = primeBuffer[i];
int exponent = exponentBuffer[i];
counts.Add(primeValue, exponent);
}
counts.Add(primeValue, exponent);
}
return true;
}
finally
{
execution.Dispose();
lease.Dispose();
}
}
finally
{
ThreadStaticPools.UlongPool.Return(primeBufferArray);
ThreadStaticPools.IntPool.Return(exponentBufferArray);
}
}
lease.Dispose();
ThreadStaticPools.UlongPool.Return(primeBufferArray);
ThreadStaticPools.IntPool.Return(exponentBufferArray);
return true;
}
}
@@ -185,55 +185,49 @@ internal static partial class PrimeOrderCalculator
private static UInt128 ExponentLoweringWide(UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
{
ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool;
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
int length = factors.Count;
// TODO: This should never trigger from production code - check
// if (length == 0)
// {
// return order;
// }
FactorEntry128[]? tempArray = pool.Rent(length + 1);
try
ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool;
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
int length = factors.Count;
// TODO: This should never trigger from production code - check
// if (length == 0)
// {
// return order;
// }
FactorEntry128[] tempArray = pool.Rent(length + 1);
Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
factorSpan.CopyTo(buffer);
if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor))
{
buffer[length] = new FactorEntry128(factors.Cofactor, 1);
length++;
}
Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
factorSpan.CopyTo(buffer);
buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value));
if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor))
for (int i = 0; i < length; i++)
{
UInt128 primeFactor = buffer[i].Value;
int exponent = buffer[i].Exponent;
for (int iteration = 0; iteration < exponent; iteration++)
{
buffer[length] = new FactorEntry128(factors.Cofactor, 1);
length++;
}
buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value));
for (int i = 0; i < length; i++)
{
UInt128 primeFactor = buffer[i].Value;
int exponent = buffer[i].Exponent;
for (int iteration = 0; iteration < exponent; iteration++)
if ((order % primeFactor) == UInt128.Zero)
{
if ((order % primeFactor) == UInt128.Zero)
UInt128 reduced = order / primeFactor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
UInt128 reduced = order / primeFactor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
order = reduced;
continue;
}
order = reduced;
continue;
}
break;
}
}
return order;
}
finally
{
pool.Return(tempArray, clearArray: false);
break;
}
}
pool.Return(tempArray, clearArray: false);
return order;
}
private static bool TryConfirmOrderWide(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
@@ -956,7 +950,7 @@ internal static partial class PrimeOrderCalculator
if (value <= ulong.MaxValue)
{
// return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None);
return PrimeTester.IsPrimeCpu((ulong)value, CancellationToken.None);
return PrimeTester.IsPrime((ulong)value);
}
if ((value & UInt128.One) == UInt128.Zero)
File diff suppressed because it is too large. Load diff
+6 -146
View File
@@ -1,151 +1,11 @@
using System;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
public static class PrimesGenerator
{
public static readonly uint[] SmallPrimes;
public static readonly ulong[] SmallPrimesPow2;
public static readonly uint[] SmallPrimesLastOne;
public static readonly ulong[] SmallPrimesPow2LastOne;
public static readonly uint[] SmallPrimesLastSeven;
public static readonly ulong[] SmallPrimesPow2LastSeven;
static PrimesGenerator()
{
BuildSmallPrimes(
out SmallPrimes,
out SmallPrimesPow2,
out SmallPrimesLastOne,
out SmallPrimesPow2LastOne,
out SmallPrimesLastSeven,
out SmallPrimesPow2LastSeven);
}
private static void BuildSmallPrimes(
out uint[] all,
out ulong[] allPow2,
out uint[] lastOne,
out ulong[] lastOnePow2,
out uint[] lastSeven,
out ulong[] lastSevenPow2)
{
uint target = PerfectNumberConstants.PrimesLimit;
int targetInt = checked((int)target);
all = new uint[targetInt];
allPow2 = new ulong[targetInt];
lastOne = new uint[targetInt];
lastOnePow2 = new ulong[targetInt];
lastSeven = new uint[targetInt];
lastSevenPow2 = new ulong[targetInt];
List<uint> primes = new(targetInt * 4 / 3 + 16);
uint candidate = 2U;
int allCount = 0;
int lastOneCount = 0;
int lastSevenCount = 0;
while (allCount < targetInt || lastOneCount < targetInt || lastSevenCount < targetInt)
{
bool isPrime = true;
int primesCount = primes.Count;
for (int i = 0; i < primesCount; i++)
{
uint p = primes[i];
// TODO: Reuse a single squared local for `p` within this loop so we avoid recalculating
// `p * p` repeatedly without storing every squared prime globally (the cache must remain
// limited to the small-divisor tables to honor the memory constraints).
if (p * p > candidate)
{
break;
}
// TODO: Replace this trial-division `%` with the sieve-based generator that avoids
// per-candidate modulo work so building the small-prime tables stops dominating
// startup time for large scans.
if (candidate % p == 0U)
{
isPrime = false;
break;
}
}
if (isPrime)
{
primes.Add(candidate);
if (allCount < targetInt)
{
// TODO: Reuse a single cached `candidateSquared` per iteration so we avoid issuing
// three separate 64-bit multiplications when populating the pow2 arrays here and
// in the last-one/last-seven branches below.
all[allCount] = candidate;
allPow2[allCount] = checked((ulong)candidate * candidate);
allCount++;
}
if (candidate != 2U)
{
LastDigit lastDigit = GetLastDigit(candidate);
if (lastOneCount < targetInt && IsAllowedForLastOne(candidate, lastDigit))
{
lastOne[lastOneCount] = candidate;
lastOnePow2[lastOneCount] = checked((ulong)candidate * candidate);
lastOneCount++;
}
if (lastSevenCount < targetInt && IsAllowedForLastSeven(candidate, lastDigit))
{
lastSeven[lastSevenCount] = candidate;
lastSevenPow2[lastSevenCount] = checked((ulong)candidate * candidate);
lastSevenCount++;
}
}
}
// TODO: Switch this increment to the Mod6 stepping table identified in the Mod6ComparisonBenchmarks
// so we skip the composite residues without relying on per-iteration `%` checks.
candidate = candidate == 2U ? 3U : candidate + 2U;
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static LastDigit GetLastDigit(uint prime)
{
return (prime % 10U) switch
{
1U => LastDigit.One,
3U => LastDigit.Three,
5U => LastDigit.Five,
7U => LastDigit.Seven,
9U => LastDigit.Nine,
_ => throw new ArgumentOutOfRangeException(nameof(prime)),
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastOne(uint prime, LastDigit lastDigit)
{
// TODO: Swap the `% 10` usage for ULongExtensions.Mod10 so the hot classification path
// reuses the benchmarked residue helper instead of repeated divisions.
return lastDigit switch
{
LastDigit.One or LastDigit.Three or LastDigit.Nine => true,
LastDigit.Seven => prime == 7U,
_ => prime == 11U,
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsAllowedForLastSeven(uint prime, LastDigit lastDigit)
{
// TODO: Route this `% 10` classification through ULongExtensions.Mod10 to match the faster
// residue helper used elsewhere in the scanner.
return lastDigit switch
{
LastDigit.Three or LastDigit.Seven or LastDigit.Nine => true,
_ => prime == 11U,
};
}
public static readonly uint[] SmallPrimes = DivisorGenerator.SmallPrimes;
public static readonly ulong[] SmallPrimesPow2 = DivisorGenerator.SmallPrimesPow2;
public static readonly uint[] SmallPrimesLastOne = DivisorGenerator.SmallPrimesLastOne;
public static readonly ulong[] SmallPrimesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne;
public static readonly uint[] SmallPrimesLastSeven = DivisorGenerator.SmallPrimesLastSeven;
public static readonly ulong[] SmallPrimesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven;
}
+150 -34
View File
@@ -6,7 +6,7 @@ using System.Threading.Tasks;
namespace PerfectNumbers.Core;
internal sealed class TaskThreadPool : IDisposable
internal sealed class TaskThreadPool
{
private readonly struct WorkItem
{
@@ -19,10 +19,15 @@ internal sealed class TaskThreadPool : IDisposable
}
private readonly ConcurrentQueue<WorkItem> _pendingTaskQueue = new();
private readonly BlockingCollection<WorkItem> _pendingTasks;
private readonly Thread[] _threads;
private readonly HashSet<Thread> _runningThreads = new();
private readonly object _lock = new();
private readonly Action<Task> _taskExecutor;
private int _disposed;
private readonly int _maximumConcurrency;
private readonly SemaphoreSlim _workAvailable = new(0);
private readonly CancellationTokenSource _disposeCancellation = new();
private int _pendingTaskCount;
private int _currentThreads;
private int _threadId;
public TaskThreadPool(int minimumThreads, Action<Task> taskExecutor)
{
@@ -32,23 +37,28 @@ internal sealed class TaskThreadPool : IDisposable
}
_taskExecutor = taskExecutor ?? throw new ArgumentNullException(nameof(taskExecutor));
_pendingTasks = new BlockingCollection<WorkItem>(_pendingTaskQueue);
_threads = new Thread[minimumThreads];
_maximumConcurrency = minimumThreads;
for (int i = 0; i < minimumThreads; i++)
List<Thread> threadsToStart = new(minimumThreads);
lock (_lock)
{
Thread thread = new(WorkerLoop)
for (int i = 0; i < minimumThreads; i++)
{
IsBackground = true,
Name = $"{nameof(TaskThreadPool)}-{i + 1}"
};
Thread thread = CreateWorkerThread();
_runningThreads.Add(thread);
_currentThreads++;
threadsToStart.Add(thread);
}
}
_threads[i] = thread;
foreach (Thread thread in threadsToStart)
{
thread.Start();
}
}
public int ThreadCount => _threads.Length;
public int ThreadCount => _maximumConcurrency;
public void Queue(Task task)
{
@@ -57,22 +67,29 @@ internal sealed class TaskThreadPool : IDisposable
throw new ArgumentNullException(nameof(task));
}
if (Volatile.Read(ref _disposed) != 0)
List<Thread>? threadsToStart;
lock (_lock)
{
throw new ObjectDisposedException(nameof(TaskThreadPool));
if (_disposeCancellation.IsCancellationRequested)
{
throw new ObjectDisposedException(nameof(TaskThreadPool));
}
_pendingTaskQueue.Enqueue(new WorkItem(task));
_pendingTaskCount++;
threadsToStart = EnsureWorkerCapacityLocked(null);
}
try
_workAvailable.Release();
if (threadsToStart is not null)
{
_pendingTasks.Add(new WorkItem(task));
}
catch (InvalidOperationException ex) when (Volatile.Read(ref _disposed) != 0)
{
throw new ObjectDisposedException(nameof(TaskThreadPool), ex);
}
catch (ObjectDisposedException ex)
{
throw new ObjectDisposedException(nameof(TaskThreadPool), ex);
foreach (Thread thread in threadsToStart)
{
thread.Start();
}
}
}
@@ -97,33 +114,132 @@ internal sealed class TaskThreadPool : IDisposable
public void Dispose()
{
if (Interlocked.Exchange(ref _disposed, 1) != 0)
Thread[] threadsToJoin;
lock (_lock)
{
return;
_disposeCancellation.Cancel();
threadsToJoin = new Thread[_runningThreads.Count];
_runningThreads.CopyTo(threadsToJoin);
}
_pendingTasks.CompleteAdding();
if (threadsToJoin.Length > 0)
{
_workAvailable.Release(threadsToJoin.Length);
}
foreach (Thread thread in _threads)
foreach (Thread thread in threadsToJoin)
{
thread.Join();
}
_pendingTasks.Dispose();
_workAvailable.Dispose();
_disposeCancellation.Dispose();
lock (_lock)
{
_runningThreads.Clear();
while (_pendingTaskQueue.TryDequeue(out _))
{
}
_pendingTaskCount = 0;
_currentThreads = 0;
}
}
private bool TryDequeueWork(out WorkItem workItem)
{
CancellationToken cancellationToken = _disposeCancellation.Token;
while (true)
{
try
{
_workAvailable.Wait(cancellationToken);
}
catch (OperationCanceledException)
{
workItem = default;
return false;
}
if (_pendingTaskQueue.TryDequeue(out workItem))
{
Interlocked.Decrement(ref _pendingTaskCount);
return true;
}
}
}
private void WorkerLoop()
{
try
{
foreach (WorkItem workItem in _pendingTasks.GetConsumingEnumerable())
if (!TryDequeueWork(out WorkItem workItem))
{
_taskExecutor(workItem.Task);
return;
}
_taskExecutor(workItem.Task);
}
catch (ObjectDisposedException)
finally
{
// The collection was disposed while the worker was waiting. Exit gracefully.
CompleteWorker();
}
}
private void CompleteWorker()
{
List<Thread>? threadsToStart = null;
lock (_lock)
{
_runningThreads.Remove(Thread.CurrentThread);
_currentThreads--;
if (!_disposeCancellation.IsCancellationRequested)
{
threadsToStart = EnsureWorkerCapacityLocked(null);
}
}
if (threadsToStart is not null)
{
foreach (Thread thread in threadsToStart)
{
thread.Start();
}
}
}
private List<Thread>? EnsureWorkerCapacityLocked(List<Thread>? threadsToStart)
{
int pendingTasks = Math.Min(_maximumConcurrency, Volatile.Read(ref _pendingTaskCount));
while (!_disposeCancellation.IsCancellationRequested && _currentThreads < pendingTasks)
{
Thread thread = CreateWorkerThread();
_runningThreads.Add(thread);
_currentThreads++;
threadsToStart ??= new List<Thread>();
threadsToStart.Add(thread);
}
return threadsToStart;
}
private Thread CreateWorkerThread()
{
int threadIndex = Interlocked.Increment(ref _threadId);
return new Thread(WorkerLoop)
{
IsBackground = true,
Name = $"{nameof(TaskThreadPool)}-{threadIndex}"
};
}
}
+1 -1
View File
@@ -6,7 +6,7 @@ namespace PerfectNumbers.Core;
/// <summary>
/// Provides thread-safe file writing for both text and binary data.
/// </summary>
public sealed class TextFileWriter : IDisposable
public sealed class TextFileWriter
{
private readonly FileStream _stream;
private StreamWriter? _textWriter;
-37
View File
@@ -260,43 +260,6 @@ namespace PerfectNumbers.Core
pool.Add(stack);
}
[ThreadStatic]
private static List<Dictionary<ulong, int>>? _ulongIntDictionaryPool;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static Dictionary<ulong, int> RentUlongIntDictionary(int capacityHint)
{
List<Dictionary<ulong, int>>? pool = _ulongIntDictionaryPool;
if (pool is not null)
{
int poolCount = pool.Count;
if (poolCount > 0)
{
int lastIndex = poolCount - 1;
Dictionary<ulong, int> dictionary = pool[lastIndex];
pool.RemoveAt(lastIndex);
dictionary.EnsureCapacity(capacityHint);
return dictionary;
}
}
return new Dictionary<ulong, int>(capacityHint);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void ReturnUlongIntDictionary(Dictionary<ulong, int> dictionary)
{
List<Dictionary<ulong, int>>? pool = _ulongIntDictionaryPool;
if (pool is null)
{
pool = new List<Dictionary<ulong, int>>(4);
_ulongIntDictionaryPool = pool;
}
pool.Add(dictionary);
}
[ThreadStatic]
private static Dictionary<ulong, int>? _factorCountDictionary;
+5 -5
View File
@@ -40,7 +40,7 @@ public static partial class ULongExtensions
{
if (((exponent >> index) & 1UL) == 0UL)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
index--;
continue;
}
@@ -59,14 +59,14 @@ public static partial class ULongExtensions
int windowLength = index - windowStart + 1;
for (int square = 0; square < windowLength; square++)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
}
ulong mask = (1UL << windowLength) - 1UL;
ulong windowValue = (exponent >> windowStart) & mask;
int tableIndex = (int)((windowValue - 1UL) >> 1);
ulong multiplier = oddPowers[tableIndex];
result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(multiplier, modulus, nPrime);
index = windowStart - 1;
}
@@ -75,7 +75,7 @@ public static partial class ULongExtensions
if (!keepMontgomery)
{
result = result.MontgomeryMultiply(1UL, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
return result;
@@ -96,7 +96,7 @@ public static partial class ULongExtensions
for (int i = 1; i < oddPowers.Length; i++)
{
previous = oddPowers[i - 1];
oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
oddPowers[i] = previous.MontgomeryMultiplyCpu(square, modulus, nPrime);
}
}
+47 -9
View File
@@ -46,6 +46,45 @@ public static partial class ULongExtensions
return state.MulMod(right, modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static ulong MontgomeryMultiplyGpu(this ulong a, ulong b, ulong modulus, ulong nPrime)
{
// Mirror the struct-based MontMul64 helper measured as the second-fastest GPU-compatible
// option in GpuUInt128Montgomery64Benchmarks so accelerator kernels reuse the validated path.
MultiplyPartsGpu(a, b, out ulong productHigh, out ulong productLow);
MultiplyPartsGpu(productLow, nPrime, out _, out ulong mLow);
MultiplyPartsGpu(mLow, modulus, out ulong mTimesModulusHigh, out ulong mTimesModulusLow);
ulong sumLow = unchecked(productLow + mTimesModulusLow);
ulong carry = sumLow < productLow ? 1UL : 0UL;
ulong sumHigh = unchecked(productHigh + mTimesModulusHigh + carry);
if (sumHigh >= modulus)
{
sumHigh -= modulus;
}
return sumHigh;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void MultiplyPartsGpu(ulong left, ulong right, out ulong high, out ulong low)
{
ulong leftLow = (uint)left;
ulong leftHigh = left >> 32;
ulong rightLow = (uint)right;
ulong rightHigh = right >> 32;
ulong lowProduct = unchecked(leftLow * rightLow);
ulong cross1 = unchecked(leftHigh * rightLow);
ulong cross2 = unchecked(leftLow * rightHigh);
ulong highProduct = unchecked(leftHigh * rightHigh);
ulong carry = unchecked((lowProduct >> 32) + (uint)cross1 + (uint)cross2);
low = unchecked((lowProduct & 0xFFFFFFFFUL) | (carry << 32));
high = unchecked(highProduct + (cross1 >> 32) + (cross2 >> 32) + (carry >> 32));
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus)
{
@@ -204,14 +243,15 @@ public static partial class ULongExtensions
ulong modulus = divisor.Modulus;
ulong nPrime = divisor.NPrime;
ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent);
return MontgomeryMultiply(montgomeryResult, 1UL, modulus, nPrime);
return MontgomeryMultiplyGpu(montgomeryResult, 1UL, modulus, nPrime);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorData divisor, ulong exponent)
{
ulong modulus = divisor.Modulus;
// Reuse the optimized MontgomeryMultiply extension highlighted in MontgomeryMultiplyBenchmarks so the GPU path matches the fastest CPU reduction.
// Reuse the GPU-compatible MontgomeryMultiply helper highlighted in MontgomeryMultiplyBenchmarks so the accelerator path
// mirrors the fastest CPU reduction without relying on UInt128 intrinsics.
// We shouldn't hit it in production code. If we do, we're doing something wrong.
// if (exponent == 0UL)
// {
@@ -233,7 +273,7 @@ public static partial class ULongExtensions
while (index >= 0)
{
ulong currentBit = (exponent >> index) & 1UL;
ulong squared = MontgomeryMultiply(result, result, modulus, nPrime);
ulong squared = MontgomeryMultiplyGpu(result, result, modulus, nPrime);
bool processWindow = currentBit != 0UL;
result = processWindow ? result : squared;
@@ -253,11 +293,11 @@ public static partial class ULongExtensions
int windowLength = processWindow ? index - windowStart + 1 : 0;
for (int square = 0; square < windowLength; square++)
{
result = MontgomeryMultiply(result, result, modulus, nPrime);
result = MontgomeryMultiplyGpu(result, result, modulus, nPrime);
}
result = processWindow
? MontgomeryMultiply(
? MontgomeryMultiplyGpu(
result,
ComputeMontgomeryOddPowerGpu(
(exponent >> windowStart) & ((1UL << windowLength) - 1UL),
@@ -330,7 +370,7 @@ public static partial class ULongExtensions
while (remaining != 0UL)
{
ulong bit = remaining & 1UL;
ulong multiplied = MontgomeryMultiply(power, baseValue, modulus, nPrime);
ulong multiplied = MontgomeryMultiplyGpu(power, baseValue, modulus, nPrime);
ulong mask = (ulong)-(long)bit;
power = (power & ~mask) | (multiplied & mask);
@@ -340,7 +380,7 @@ public static partial class ULongExtensions
break;
}
baseValue = MontgomeryMultiply(baseValue, baseValue, modulus, nPrime);
baseValue = MontgomeryMultiplyGpu(baseValue, baseValue, modulus, nPrime);
}
return power;
@@ -417,7 +457,6 @@ public static partial class ULongExtensions
ulong result = 0UL;
GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
// Keep this commented out. It should never happen in production code.
@@ -444,7 +483,6 @@ public static partial class ULongExtensions
resultBuffer.Dispose();
exponentBuffer.Dispose();
execution.Dispose();
lease.Dispose();
// Keep this commented. We don't want to catch any exceptions. All should crash the scanner.
// catch (AcceleratorException)
+4 -4
View File
@@ -553,7 +553,7 @@ public static partial class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)
public static ulong MontgomeryMultiplyCpu(this ulong a, ulong b, ulong modulus, ulong nPrime)
{
ulong tLow = unchecked(a * b);
ulong m = unchecked(tLow * nPrime);
@@ -581,7 +581,7 @@ public static partial class ULongExtensions
{
if ((remainingExponent & 1UL) != 0UL)
{
result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
}
remainingExponent >>= 1;
@@ -590,10 +590,10 @@ public static partial class ULongExtensions
break;
}
baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
baseVal = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
}
return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime);
return keepMontgomery ? result : result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
}
+2 -2
View File
@@ -131,12 +131,12 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels,
* [done] The CPU heuristic now operates as a pure trial-division sweep without maintaining per-call summaries, relying on `PrepareHeuristicDivisor` only when downstream consumers request Montgomery data.
* [done] Introduced the temporary `UseHeuristicGroupBTrialDivision` gate so current builds execute Group A locally and then fall back to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, keeping the full Group B implementation available for reactivation.
7. [done] **Implement `IsPrimeGpu`**
* [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `sqrtLimit`/`nMod10` inputs, and falling back to CPU when GPUs are disabled.
* [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `maxDivisorSquare`/`nMod10` inputs, and falling back to CPU when GPUs are disabled.
* [done] Added `PrimeTesterKernels.HeuristicTrialDivisionKernel` so GPU batches respect Group A/B ordering while emitting hit flags for early exits.
* [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic.
* [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement.
8. [done] **Refactor by-divisor CPU/GPU testers to reuse heuristics**
* [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop while preserving the existing status counters.
* [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop and retiring the legacy status counter.
* [done] Streamed GPU batches from the same enumerator and Montgomery preparation data, removing the filtered-divisor scratch arrays while continuing to feed the existing kernel and hit bookkeeping.
9. [done] **Route legacy prime checks through heuristics**
* [done] Introduce `HeuristicPrimeTester` so heuristic CPU/GPU checks live alongside the legacy `PrimeTester` implementations while callers opt in explicitly.
+8
View File
@@ -0,0 +1,8 @@
# GPU context pooling and heuristic divisor updates
- [x] Review existing GPU context pooling in `PerfectNumbers.Core/PrimeTester.cs` and determine shared table allocations that can be scoped per trailing digit.
- [x] Split `PrimeTesterGpuContextPool` into trailing-digit specific pools so each lease only uploads the divisor tables required for that digit.
- [x] Update GPU warm-up and rental code paths to pick pools based on `n % 10` and ensure thread-safe reuse.
- [x] Replace square-table usage in GPU heuristic trial division with square-root limits so divisor tables no longer need squared values on the accelerator.
- [x] Remove GPU square buffers and adjust kernels, table types, and callers accordingly.
- [x] Build solution to confirm the new pooling and divisor logic compiles; tests documented in summary.
-1
View File
@@ -41,7 +41,6 @@ Montgomery Primitives
Cycle-Aware Modular Arithmetic
Status & Progress Reporting
- PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs:296: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
Cycle Remainder Steppers
- PerfectNumbers.Core/CycleRemainderStepper.cs:28: 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.