Compare commits

..
Author SHA1 Message Date
Marek Lemańczyk 0466d98434 Hoist by-divisor invariants and split mod-8 updates 2025-10-27 16:48:02 +01:00
64 changed files with 2951 additions and 4989 deletions

No files matched your search

@@ -32,7 +32,7 @@ public class GpuUInt128Montgomery64Benchmarks
[Benchmark(Baseline = true)]
public ulong ExtensionMontgomeryMultiply()
{
return Input.Left.MontgomeryMultiplyCpu(Input.Right, Input.Modulus, Input.NPrime);
return Input.Left.MontgomeryMultiply(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.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
@@ -736,22 +736,22 @@ public class Pow2ModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent);
for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--)
{
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(result, modulus, nPrime);
if (((exponent >> bitIndex) & 1UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
}
}
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime);
ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime);
if ((nibble & 1UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime);
result = result.MontgomeryMultiply(base0, modulus, nPrime);
}
if ((nibble & 2UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime);
result = result.MontgomeryMultiply(base1, modulus, nPrime);
}
if ((nibble & 4UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime);
result = result.MontgomeryMultiply(base2, modulus, nPrime);
}
if ((nibble & 8UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base3, modulus, nPrime);
result = result.MontgomeryMultiply(base3, modulus, nPrime);
}
baseVal = base3.MontgomeryMultiplyCpu(base3, modulus, nPrime);
baseVal = base3.MontgomeryMultiply(base3, modulus, nPrime);
remainingExponent >>= 4;
}
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor)
@@ -817,16 +817,16 @@ public class Pow2ModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
const int WindowSize = 5;
Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)];
oddPowers[0] = baseVal;
ulong baseSquared = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
ulong baseSquared = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
for (int i = 1; i < oddPowers.Length; i++)
{
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiplyCpu(baseSquared, modulus, nPrime);
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(baseSquared, modulus, nPrime);
}
int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent);
@@ -835,7 +835,7 @@ public class Pow2ModBenchmarks
{
if (((exponent >> bitIndex) & 1UL) == 0UL)
{
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(result, modulus, nPrime);
}
result = result.MontgomeryMultiplyCpu(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
result = result.MontgomeryMultiply(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
bitIndex = windowStart - 1;
}
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
@@ -300,22 +300,22 @@ public class Pow2MontgomeryModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent);
for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--)
{
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(result, modulus, nPrime);
if (((exponent >> bitIndex) & 1UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
}
}
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime);
ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime);
if ((nibble & 1UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime);
result = result.MontgomeryMultiply(base0, modulus, nPrime);
}
if ((nibble & 2UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime);
result = result.MontgomeryMultiply(base1, modulus, nPrime);
}
if ((nibble & 4UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime);
result = result.MontgomeryMultiply(base2, modulus, nPrime);
}
if ((nibble & 8UL) != 0UL)
{
result = result.MontgomeryMultiplyCpu(base3, modulus, nPrime);
result = result.MontgomeryMultiply(base3, modulus, nPrime);
}
baseVal = base3.MontgomeryMultiplyCpu(base3, modulus, nPrime);
baseVal = base3.MontgomeryMultiply(base3, modulus, nPrime);
remainingExponent >>= 4;
}
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
// TODO: Update it with the latest implementations
@@ -382,16 +382,16 @@ public class Pow2MontgomeryModBenchmarks
if (exponent == 0UL)
{
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
const int WindowSize = 5;
Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)];
oddPowers[0] = baseVal;
ulong baseSquared = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
ulong baseSquared = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
for (int i = 1; i < oddPowers.Length; i++)
{
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiplyCpu(baseSquared, modulus, nPrime);
oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(baseSquared, modulus, nPrime);
}
int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent);
@@ -400,7 +400,7 @@ public class Pow2MontgomeryModBenchmarks
{
if (((exponent >> bitIndex) & 1UL) == 0UL)
{
result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(result, modulus, nPrime);
}
result = result.MontgomeryMultiplyCpu(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
result = result.MontgomeryMultiply(oddPowers[(int)(windowValue >> 1)], modulus, nPrime);
bitIndex = windowStart - 1;
}
return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return result.MontgomeryMultiply(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 static readonly HeuristicCombinedPrimeTester _heuristicCombinedTester = new();
private static readonly HeuristicPrimeTester _heuristicTester = new();
private HeuristicPrimeTester _cpuTester = null!;
private HeuristicPrimeTester _gpuTester = null!;
private ulong[] _candidates = [];
[ParamsSource(nameof(GetBenchmarkCases))]
@@ -22,6 +22,8 @@ public class PrimeTesterBenchmarks
[GlobalSetup]
public void GlobalSetup()
{
_cpuTester = new HeuristicPrimeTester();
_gpuTester = new HeuristicPrimeTester();
_candidates = BenchmarkCase.Candidates;
if (_candidates.Length == 0)
@@ -29,21 +31,23 @@ public class PrimeTesterBenchmarks
return;
}
_ = _heuristicTester.IsPrimeCpu(_candidates[0]);
_ = _heuristicTester.IsPrimeGpu(_candidates[0]);
CancellationToken cancellationToken = CancellationToken.None;
_cpuTester.IsPrimeCpu(_candidates[0], cancellationToken);
_gpuTester.IsPrimeGpu(_candidates[0]);
_ = Prime.Numbers.IsPrime(_candidates[0]);
}
[Benchmark(Baseline = true)]
public int HeuristicCpu()
{
HeuristicPrimeTester tester = _heuristicTester;
HeuristicPrimeTester tester = _cpuTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (tester.IsPrimeCpu(values[i]))
if (tester.IsPrimeCpu(values[i], cancellationToken))
{
primeCount++;
}
@@ -55,8 +59,9 @@ public class PrimeTesterBenchmarks
[Benchmark]
public int HeuristicGpu()
{
HeuristicPrimeTester tester = _heuristicTester;
HeuristicPrimeTester tester = _gpuTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
@@ -70,52 +75,18 @@ public class PrimeTesterBenchmarks
return primeCount;
}
[Benchmark]
public int HeuristicCombinedGpu()
{
ulong[] values = _candidates;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
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]
// This benchmark is very slow. No point in comparing it.
// [Benchmark]
public int NonHeuristicGpu()
{
PrimeTester tester = _nonHeuristicTester;
{
// PrimeTester tester = _gpuTester;
ulong[] values = _candidates;
CancellationToken cancellationToken = CancellationToken.None;
int primeCount = 0;
for (int i = 0; i < values.Length; i++)
{
if (tester.IsPrimeGpu(values[i]))
if (PrimeTester.Exclusive.IsPrimeGpu(values[i], cancellationToken))
{
primeCount++;
}
@@ -127,12 +98,14 @@ 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.IsPrime(values[i]))
if (PrimeTester.IsPrimeCpu(values[i], cancellationToken))
{
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)
{
PrimeTester.IsPrime(n).Should().Be(expected);
new PrimeTester().IsPrime(n, CancellationToken.None).Should().Be(expected);
}
[Fact]
+34 -45
View File
@@ -202,15 +202,14 @@ internal readonly struct CliArguments
if (argument.StartsWith("--prime=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> primeValue = argument.AsSpan("--prime=".Length);
startPrime = Utf8CliParser.ParseUInt64(primeValue);
startPrime = Utf8CliParser.ParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1));
startPrimeProvided = true;
continue;
}
if (argument.StartsWith("--max-prime=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan("--max-prime=".Length);
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
if (UInt128.TryParse(value, NumberStyles.None, CultureInfo.InvariantCulture, out UInt128 parsedMaxPrime))
{
maxPrimeLimit = parsedMaxPrime;
@@ -239,13 +238,13 @@ internal readonly struct CliArguments
if (argument.StartsWith("--threads=", StringComparison.OrdinalIgnoreCase))
{
threadCount = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--threads=".Length)));
threadCount = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan("--mersenne=".Length);
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
if (value.Equals("pow2mod", StringComparison.OrdinalIgnoreCase))
{
kernelType = GpuKernelType.Pow2Mod;
@@ -282,7 +281,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--divisor-cycles-limit=".Length), out ulong parsedLimit))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedLimit))
{
divisorCyclesSearchLimit = parsedLimit;
}
@@ -292,7 +291,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--residue-max-k=".Length), out ulong parsedResidueMax))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedResidueMax))
{
residueKMax = parsedResidueMax;
}
@@ -302,7 +301,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase))
{
useMersenneOnGpu = !argument.AsSpan("--mersenne-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useMersenneOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
@@ -320,7 +319,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--order-warmup-limit=".Length), out ulong parsedLimit))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedLimit))
{
orderWarmupLimitOverride = parsedLimit;
}
@@ -336,14 +335,14 @@ internal readonly struct CliArguments
if (argument.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan("--ntt=".Length);
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
nttBackend = value.Equals("staged", StringComparison.OrdinalIgnoreCase) ? NttBackend.Staged : NttBackend.Reference;
continue;
}
if (argument.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan("--mod-reduction=".Length);
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase))
{
modReductionMode = ModReductionMode.Mont64;
@@ -366,25 +365,25 @@ internal readonly struct CliArguments
if (argument.StartsWith("--primes-device=", StringComparison.OrdinalIgnoreCase))
{
forcePrimeKernelsOnCpu = argument.AsSpan("--primes-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
forcePrimeKernelsOnCpu = argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--order-device=", StringComparison.OrdinalIgnoreCase))
{
useOrderOnGpu = !argument.AsSpan("--order-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useOrderOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--rle-blacklist=", StringComparison.OrdinalIgnoreCase))
{
rleBlacklistPath = argument["--rle-blacklist=".Length..];
rleBlacklistPath = argument[(argument.IndexOf('=') + 1)..];
continue;
}
if (argument.StartsWith("--rle-hard-max=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--rle-hard-max=".Length), out ulong parsedMax))
if (Utf8CliParser.TryParseUInt64(argument.AsSpan(argument.IndexOf('=') + 1), out ulong parsedMax))
{
rleHardMaxP = parsedMax;
}
@@ -394,14 +393,14 @@ internal readonly struct CliArguments
if (argument.StartsWith("--rle-only-last7=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan("--rle-only-last7=".Length);
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
rleOnlyLast7 = !value.Equals("false", StringComparison.OrdinalIgnoreCase) && !value.Equals("0", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--zero-hard=", StringComparison.OrdinalIgnoreCase))
{
if (Utf8CliParser.TryParseDouble(argument.AsSpan("--zero-hard=".Length), out double zeroFraction))
if (Utf8CliParser.TryParseDouble(argument.AsSpan(argument.IndexOf('=') + 1), out double zeroFraction))
{
zeroFractionHard = zeroFraction;
}
@@ -411,7 +410,7 @@ internal readonly struct CliArguments
if (argument.StartsWith("--zero-conj=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = argument.AsSpan("--zero-conj=".Length);
ReadOnlySpan<char> value = argument.AsSpan(argument.IndexOf('=') + 1);
int colon = value.IndexOf(':');
if (colon > 0)
{
@@ -428,97 +427,87 @@ internal readonly struct CliArguments
if (argument.StartsWith("--lucas=", StringComparison.OrdinalIgnoreCase))
{
useMersenneOnGpu = !argument.AsSpan("--lucas=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useMersenneOnGpu = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase))
if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase) ||
argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase) ||
argument.StartsWith("--accelerator=", 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);
forcePrimeKernelsOnCpu = argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--results-dir=", StringComparison.OrdinalIgnoreCase))
{
resultsDirectory = argument["--results-dir=".Length..];
resultsDirectory = argument[(argument.IndexOf('=') + 1)..];
continue;
}
if (argument.StartsWith("--results-prefix=", StringComparison.OrdinalIgnoreCase))
{
resultsPrefix = argument["--results-prefix=".Length..];
resultsPrefix = argument[(argument.IndexOf('=') + 1)..];
continue;
}
if (argument.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase))
{
gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-threads=".Length)));
gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--gpu-prime-batch=", StringComparison.OrdinalIgnoreCase))
{
gpuPrimeBatch = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-batch=".Length)));
gpuPrimeBatch = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--ll-slice=", StringComparison.OrdinalIgnoreCase))
{
sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--ll-slice=".Length)));
sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--gpu-scan-batch=", StringComparison.OrdinalIgnoreCase))
{
scanBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-scan-batch=".Length)));
scanBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--block-size=", StringComparison.OrdinalIgnoreCase))
{
blockSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--block-size=".Length)));
blockSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--filter-p=", StringComparison.OrdinalIgnoreCase))
{
filterFile = argument["--filter-p=".Length..];
filterFile = argument[(argument.IndexOf('=') + 1)..];
continue;
}
if (argument.StartsWith("--write-batch-size=", StringComparison.OrdinalIgnoreCase))
{
writeBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--write-batch-size=".Length)));
writeBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
continue;
}
if (argument.StartsWith("--divisor-cycles=", StringComparison.OrdinalIgnoreCase))
{
cyclesPath = argument["--divisor-cycles=".Length..];
cyclesPath = argument[(argument.IndexOf('=') + 1)..];
continue;
}
if (argument.StartsWith("--divisor-cycles-device=", StringComparison.OrdinalIgnoreCase))
{
useGpuCycles = !argument.AsSpan("--divisor-cycles-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase);
useGpuCycles = !argument.AsSpan(argument.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
continue;
}
if (argument.StartsWith("--divisor-cycles-batch=", StringComparison.OrdinalIgnoreCase))
{
cyclesBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--divisor-cycles-batch=".Length)));
cyclesBatchSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan(argument.IndexOf('=') + 1)));
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> trimmed;
ReadOnlySpan<char> span;
while ((line = reader.ReadLine()) is not null)
{
if (string.IsNullOrWhiteSpace(line))
@@ -22,67 +22,36 @@ internal static class CalculationResultsFile
continue;
}
trimmed = line.AsSpan().Trim();
if (trimmed.IsEmpty || trimmed[0] == '#')
span = line.AsSpan().Trim();
if (span.IsEmpty || span[0] == '#')
{
continue;
}
int length = trimmed.Length;
int index = 0;
if (!maxPrimeConfigured)
while (index < span.Length)
{
while (index < length)
if (!char.IsDigit(span[index]))
{
char current = trimmed[index];
uint digit = (uint)(current - '0');
if (digit > 9U)
{
index++;
continue;
}
index++;
continue;
}
int start = index++;
while (index < length && (uint)(trimmed[index] - '0') <= 9U)
{
index++;
}
int start = index;
index++;
while (index < span.Length && char.IsDigit(span[index]))
{
index++;
}
if (Utf8CliParser.TryParseUInt64(trimmed[start..index], out ulong parsed))
if (Utf8CliParser.TryParseUInt64(span[start..index], out ulong parsed))
{
if (!maxPrimeConfigured || (UInt128)parsed <= maxPrimeLimit)
{
candidates.Add(parsed);
}
}
}
else
{
UInt128 limit = maxPrimeLimit;
while (index < length)
{
char current = trimmed[index];
uint digit = (uint)(current - '0');
if (digit > 9U)
else
{
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++;
}
skippedByLimit++;
}
}
}
+19 -14
View File
@@ -28,8 +28,6 @@ internal static class Program
[ThreadStatic]
private static bool _lastCompositeP;
private static bool _runPrimesOnCpu;
private static void Main(string[] args)
{
try
@@ -51,21 +49,18 @@ internal static class Program
NttGpuMath.GpuTransformBackend = _cliArguments.NttBackend;
NttGpuMath.ReductionMode = _cliArguments.ModReductionMode;
_runPrimesOnCpu = _cliArguments.ForcePrimeKernelsOnCpu;
GpuContextPool.ForceCpu = _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;
@@ -109,7 +104,14 @@ internal static class Program
testPrimeCandidateLimit = Math.Max(1, threadCount) * 3;
}
_byDivisorStartPrime = useByDivisor && startPrimeProvided ? currentP : 0UL;
if (useByDivisor)
{
_byDivisorStartPrime = startPrimeProvided ? currentP : 0UL;
}
else
{
_byDivisorStartPrime = 0UL;
}
if (useByDivisor && !testMode && string.IsNullOrEmpty(filterFile))
{
@@ -119,7 +121,7 @@ internal static class Program
// Apply GPU prime sieve runtime configuration
GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads);
PrimeTester.GpuBatchSize = gpuPrimeBatch;
PrimeTester.GpuBatchSize = Math.Max(1, gpuPrimeBatch);
MersenneDivisorCycles mersenneDivisorCycles = new();
if (!File.Exists(cyclesPath) || continueCyclesGeneration)
@@ -175,7 +177,10 @@ internal static class Program
Console.WriteLine("Divisor cycles are ready");
blockSize = useByDivisor ? 1 : blockSize;
if (useByDivisor)
{
blockSize = 1;
}
CandidatesCalculator.Initialize(currentP);
PrimeTesters = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
@@ -261,7 +266,7 @@ internal static class Program
orderWarmupLimitOverride ?? 5_000_000UL,
NttGpuMath.ReductionMode,
mersenneOnGpu ? "gpu" : "cpu",
_runPrimesOnCpu ? "cpu" : "gpu",
GpuContextPool.ForceCpu ? "cpu" : "gpu",
orderOnGpu ? "gpu" : "cpu");
if (!string.IsNullOrEmpty(_cliArguments.ResultsPrefix))
@@ -626,9 +631,9 @@ internal static class Program
// If primes-device=gpu, route p primality through GPU-assisted sieve with deterministic MR validation.
if (!_cliArguments.UseByDivisor)
{
if (!(_runPrimesOnCpu
? PrimeTester.IsPrime(p)
: PrimeTesters.Value!.IsPrimeGpu(p)))
if (!(GpuContextPool.ForceCpu
? PrimeTesters.Value!.IsPrime(p, CancellationToken.None)
: PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None)))
{
_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<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> PrimeTesterSmallPrimeKernel =
private static readonly Action<Index1D, ArrayView<ulong>, ArrayView<uint>, 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<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterSmallPrimeKernel);
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>>(PrimeTesterSmallPrimeKernel);
_ = kernel;
}
@@ -18,33 +18,22 @@ public class GpuPrimeLimiterTests
var startedSecond = new TaskCompletionSource<bool>();
var sw = new Stopwatch();
var first = GpuPrimeWorkLimiter.Acquire();
using var first = GpuPrimeWorkLimiter.Acquire();
sw.Start();
var secondTask = Task.Run(() =>
{
var second = GpuPrimeWorkLimiter.Acquire();
try
{
startedSecond.TrySetResult(true);
}
finally
{
second.Dispose();
}
using var second = GpuPrimeWorkLimiter.Acquire();
startedSecond.TrySetResult(true);
});
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();
}
// 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();
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().IsPrime(p, true).Should().Be(expected);
new MersenneNumberLucasLehmerGpuTester().IsMersennePrime(p).Should().Be(expected);
}
[Fact]
@@ -26,7 +26,7 @@ public class LucasLehmerPrimeTesterTests
const ulong exponent = 128UL;
var tester = new MersenneNumberLucasLehmerGpuTester();
tester.WarmUpNttParameters(exponent);
tester.IsPrime(exponent, true).Should().BeFalse();
tester.IsMersennePrime(exponent).Should().BeFalse();
}
[Fact]
@@ -49,7 +49,7 @@ public class LucasLehmerPrimeTesterTests
{
var tester = new MersenneNumberLucasLehmerGpuTester();
tester.WarmUpNttParameters(exponent);
tester.IsPrime(exponent, true).Should().BeFalse();
tester.IsMersennePrime(exponent).Should().BeFalse();
}
[Theory]
@@ -10,17 +10,19 @@ 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).Should().BeTrue();
tester.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue();
}
}
finally
{
GpuContextPool.ForceCpu = false;
GpuContextPool.DisposeAll();
}
}
@@ -28,17 +30,19 @@ 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).Should().BeFalse();
tester.IsPrimeGpu(composite, CancellationToken.None).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,31 +17,19 @@ public class GpuWorkLimiterTests
var startedSecond = new TaskCompletionSource<bool>();
var sw = new Stopwatch();
var first = GpuWorkLimiter.Acquire();
using var first = GpuWorkLimiter.Acquire();
sw.Start();
var secondTask = Task.Run(() =>
{
var second = GpuWorkLimiter.Acquire();
try
{
startedSecond.TrySetResult(true);
}
finally
{
second.Dispose();
}
using var second = GpuWorkLimiter.Acquire();
startedSecond.TrySetResult(true);
});
try
{
await Task.Delay(50);
startedSecond.Task.IsCompleted.Should().BeFalse();
}
finally
{
first.Dispose();
}
await Task.Delay(50);
startedSecond.Task.IsCompleted.Should().BeFalse();
first.Dispose();
await secondTask.WaitAsync(TimeSpan.FromSeconds(1));
sw.Stop();
@@ -1,239 +0,0 @@
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,25 +30,18 @@ public class MersenneNumberDivisorCpuTesterTests
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(13UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
byte[] hits = new byte[primes.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
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 });
}
[Fact]
@@ -58,22 +51,15 @@ public class MersenneNumberDivisorCpuTesterTests
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(43UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
[Fact]
@@ -83,22 +69,15 @@ public class MersenneNumberDivisorCpuTesterTests
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(61UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
}
@@ -149,20 +149,13 @@ public class MersenneNumberDivisorGpuTesterTests
var tester = new MersenneNumberDivisorByDivisorGpuTester();
tester.ConfigureFromMaxPrime(43UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL };
byte[] hits = new byte[primes.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().ContainInOrder(new byte[] { 0, 1, 0 });
}
[Fact]
@@ -172,25 +165,18 @@ public class MersenneNumberDivisorGpuTesterTests
var tester = new MersenneNumberDivisorByDivisorGpuTester();
tester.ConfigureFromMaxPrime(43UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL, 43UL };
byte[] hits = new byte[primes.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
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 });
}
[Fact]
@@ -204,21 +190,14 @@ public class MersenneNumberDivisorGpuTesterTests
tester.ConfigureFromMaxPrime(47UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
}
[Fact]
@@ -231,23 +210,16 @@ public class MersenneNumberDivisorGpuTesterTests
};
tester.ConfigureFromMaxPrime(43UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
[Fact]
@@ -260,23 +232,16 @@ public class MersenneNumberDivisorGpuTesterTests
};
tester.ConfigureFromMaxPrime(61UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
[Fact]
@@ -332,21 +297,14 @@ public class MersenneNumberDivisorGpuTesterTests
tester.ConfigureFromMaxPrime(47UL);
var session = tester.CreateDivisorSession();
try
{
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
using var session = tester.CreateDivisorSession();
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 });
}
finally
{
session.Dispose();
}
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
}
}
@@ -17,29 +17,12 @@ 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();
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);
tester.IsPrime(n, CancellationToken.None).Should().Be(expected);
}
[Fact]
@@ -47,14 +30,20 @@ public class PrimeTesterIsPrimeTests
public void IsPrimeGpu_falls_back_to_cpu_when_forced()
{
var tester = new PrimeTester();
PrimeTester.IsPrime(137UL).Should().BeTrue();
PrimeTester.IsPrime(341UL).Should().BeFalse();
GpuContextPool.ForceCpu = true;
tester.IsPrimeGpu(137UL, CancellationToken.None).Should().BeTrue();
tester.IsPrimeGpu(341UL, CancellationToken.None).Should().BeFalse();
GpuContextPool.ForceCpu = false;
}
[Fact]
[Trait("Category", "Fast")]
public void IsPrimeGpu_handles_parallel_calls_for_known_primes()
{
GpuContextPool.ForceCpu = false;
GpuPrimeWorkLimiter.SetLimit(64);
PrimeTester.GpuBatchSize = 256_000;
@@ -64,7 +53,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]);
results[i] = testers.Value!.IsPrimeGpu(primes[i], CancellationToken.None);
});
results.Should().AllBeEquivalentTo(true);
@@ -73,7 +62,7 @@ public class PrimeTesterIsPrimeTests
{
foreach (ulong prime in primes)
{
testers.Value!.IsPrimeGpu(prime).Should().BeTrue();
testers.Value!.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue();
}
});
}
@@ -32,39 +32,14 @@ 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++)
{
@@ -83,64 +58,6 @@ 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,8 +2,11 @@ namespace PerfectNumbers.Core.Cpu;
internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{
private bool _disposed;
public void Reset()
{
_disposed = false;
}
public void CheckDivisor(
@@ -64,6 +67,12 @@ 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);
}
}
@@ -135,78 +135,40 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
ulong stepHigh = step.High;
ulong stepLow = step.Low;
ulong limitHigh = limit.High;
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 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;
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)
{
bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0;
if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
ulong candidate = divisor.Low;
bool admissible = ((decimalMask >> remainder10) & 1) != 0
&& (remainder8 == 1 || remainder8 == 7)
&& remainder3 != 0
&& remainder5 != 0
&& remainder7 != 0
&& remainder11 != 0;
if (admissible)
{
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.
@@ -251,139 +213,18 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
divisor.Add(step);
remainder10 = AddMod10(remainder10, step10);
remainder10 = AddMod(remainder10, step10, (byte)10);
remainder8 = AddMod8(remainder8, step8);
remainder3 = AddMod3(remainder3, step3);
remainder7 = AddMod7(remainder7, step7);
remainder11 = AddMod11(remainder11, step11);
remainder5 = AddMod(remainder5, step5, (byte)5);
remainder3 = AddMod(remainder3, step3, (byte)3);
remainder7 = AddMod(remainder7, step7, (byte)7);
remainder11 = AddMod(remainder11, step11, (byte)11);
}
processedAll = true;
processedAll = divisor.CompareTo(limit) > 0;
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;
}
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)
{
ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
@@ -391,57 +232,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod3(byte value, byte delta)
private static byte AddMod(byte value, byte delta, byte modulus)
{
const int Modulus = 3;
int sum = value + delta;
if (sum >= Modulus)
if (sum >= modulus)
{
sum -= Modulus;
}
return (byte)sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod7(byte value, byte delta)
{
const int Modulus = 7;
int sum = value + delta;
if (sum >= Modulus)
{
sum -= Modulus;
}
return (byte)sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod10(byte value, byte delta)
{
const int Modulus = 10;
int sum = value + delta;
if (sum >= Modulus)
{
sum -= Modulus;
}
return (byte)sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod11(byte value, byte delta)
{
const int Modulus = 11;
int sum = value + delta;
if (sum >= Modulus)
{
sum -= Modulus;
sum -= modulus;
}
return (byte)sum;
+2
View File
@@ -227,6 +227,7 @@ 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);
@@ -308,6 +309,7 @@ public sealed class DivisorCycleCache
ArrayPool<byte>.Shared.Return(rentedStatus!, clearArray: false);
}
execution.Dispose();
gpuLease.Dispose();
}
+19 -253
View File
@@ -1,277 +1,43 @@
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);
static DivisorGenerator()
{
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
{
LastDigit.Seven => DecimalMaskWhenLastIsSeven,
LastDigit.Three => DecimalMaskWhenLastIsThree,
LastDigit.Nine => DecimalMaskWhenLastIsNine,
_ => DecimalMaskOtherwise,
};
return lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : 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)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsValidDivisor(
byte remainder10,
byte remainder8,
byte remainder3,
byte remainder5,
byte remainder7,
byte remainder11,
ushort decimalMask)
{
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)
if (((decimalMask >> remainder10) & 1) == 0)
{
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;
return false;
}
lastOneWithoutThree = lastOneWithoutThreeList.ToArray();
lastOneWithoutThreePow2 = lastOneWithoutThreePow2List.ToArray();
lastSevenWithoutThree = lastSevenWithoutThreeList.ToArray();
lastSevenWithoutThreePow2 = lastSevenWithoutThreePow2List.ToArray();
lastNineWithoutThree = lastNineWithoutThreeList.ToArray();
lastNineWithoutThreePow2 = lastNineWithoutThreePow2List.ToArray();
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static LastDigit GetLastDigit(uint prime)
{
return (prime % 10U) switch
if (remainder8 != 1 && remainder8 != 7)
{
1U => LastDigit.One,
3U => LastDigit.Three,
5U => LastDigit.Five,
7U => LastDigit.Seven,
9U => LastDigit.Nine,
_ => throw new ArgumentOutOfRangeException(nameof(prime)),
};
}
return false;
}
[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
if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0)
{
LastDigit.One or LastDigit.Three or LastDigit.Nine => true,
LastDigit.Seven => prime == 7U,
_ => prime == 11U,
};
}
return false;
}
[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,
};
return true;
}
}
@@ -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.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(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.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(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.MontgomeryMultiplyCpu(montgomeryDelta, _modulus, _nPrime);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime);
PreviousExponent = exponent;
isUnity = _currentMontgomery == _montgomeryOne;
return true;
@@ -138,5 +138,5 @@ internal struct ExponentRemainderStepperCpu
PreviousExponent = exponent;
}
private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime);
private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);
}
@@ -1,142 +0,0 @@
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();
}
}
+99 -53
View File
@@ -1,4 +1,3 @@
using System;
using System.Collections.Concurrent;
using ILGPU;
using ILGPU.Runtime;
@@ -7,19 +6,26 @@ namespace PerfectNumbers.Core.Gpu;
public static class GpuContextPool
{
// 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; }
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;
public PooledContext()
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)
{
Context = Context.CreateDefault();
Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context);
DivisorGeneratorGpuCache.WarmUpByDivisorTables(Accelerator);
Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
// 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
@@ -28,90 +34,130 @@ 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 int WarmedGpuContextCount;
private static readonly object CreationLock = new();
// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
public static GpuContextLease Rent()
{
if (GpuPool.TryDequeue(out var gpu))
{
return new GpuContextLease(gpu);
}
// Create a new accelerator when the pool does not have one available.
return new GpuContextLease(new PooledContext());
return RentPreferred(ForceCpu);
}
public static void WarmUpPool(int threadCount)
// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
public static GpuContextLease RentPreferred(bool preferCpu)
{
if (threadCount <= WarmedGpuContextCount)
if (PoolingEnabled)
{
return;
if (preferCpu && CpuPool.TryDequeue(out var cpu))
{
return new GpuContextLease(cpu, CreationLock);
}
if (!preferCpu && GpuPool.TryDequeue(out var gpu))
{
return new GpuContextLease(gpu, CreationLock);
}
}
int toCreate = threadCount - WarmedGpuContextCount;
var contexts = new PooledContext[toCreate];
for (int i = 0; i < toCreate; i++)
// Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
lock (CreationLock)
{
contexts[i] = new PooledContext();
return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
}
for (int i = 0; i < toCreate; i++)
{
GpuPool.Enqueue(contexts[i]);
}
WarmedGpuContextCount = threadCount;
}
public static void DisposeAll()
{
if (!PoolingEnabled)
{
return;
}
while (CpuPool.TryDequeue(out var cpu))
{
cpu.Dispose();
}
while (GpuPool.TryDequeue(out var gpu))
{
gpu.Dispose();
}
DivisorGeneratorGpuCache.DisposeAll();
PrimeTester.DisposeGpuContexts();
WarmedGpuContextCount = 0;
}
private static void Return(PooledContext ctx)
{
ctx.Accelerator.Synchronize();
GpuPool.Enqueue(ctx);
}
public readonly struct GpuContextLease
{
private readonly PooledContext _ctx;
internal GpuContextLease(PooledContext ctx)
if (PoolingEnabled)
{
_ctx = ctx;
ctx.Accelerator.Synchronize();
if (ctx.IsCpu)
{
CpuPool.Enqueue(ctx);
}
else
{
GpuPool.Enqueue(ctx);
}
}
public readonly Context Context => _ctx.Context;
public readonly Accelerator Accelerator => _ctx.Accelerator;
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 void Dispose()
{
Return(_ctx);
}
}
disposed = true;
}
}
}
+105 -45
View File
@@ -1,4 +1,3 @@
using System;
using System.Collections.Concurrent;
using ILGPU;
using ILGPU.Runtime;
@@ -10,26 +9,30 @@ public sealed class GpuKernelLease
{
private static readonly ConcurrentQueue<GpuKernelLease> Pool = new();
private Action? _releaseLimiter;
private IDisposable? _limiter;
private GpuContextLease? _gpu;
private KernelContainer? _kernels;
private AcceleratorStream? _stream;
private object? _executionLock;
private bool _disposed;
private GpuKernelLease()
{
}
internal static GpuKernelLease Rent(GpuPrimeWorkLimiter.Lease limiter, GpuContextLease gpu, KernelContainer kernels)
internal static GpuKernelLease Rent(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels)
{
if (!Pool.TryDequeue(out var lease))
{
lease = new GpuKernelLease();
}
lease._releaseLimiter = limiter.Dispose;
lease._limiter = limiter;
lease._gpu = gpu;
lease._kernels = kernels;
lease._stream = null;
lease._executionLock = gpu.ExecutionLock;
lease._disposed = false;
return lease;
}
@@ -51,17 +54,22 @@ 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
return KernelContainer.InitOnce(ref _kernels!.Order, () =>
lock (_gpu!.Value.KernelInitLock)
{
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>>>();
});
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>>>();
});
}
}
}
@@ -70,12 +78,15 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
return KernelContainer.InitOnce(ref _kernels!.Pow2Mod, () =>
lock (_gpu!.Value.KernelInitLock)
{
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>>>();
});
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>>>();
});
}
}
}
@@ -84,12 +95,15 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
return KernelContainer.InitOnce(ref _kernels!.Incremental, () =>
lock (_gpu!.Value.KernelInitLock)
{
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>>>();
});
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>>>();
});
}
}
}
@@ -98,12 +112,15 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
return KernelContainer.InitOnce(ref _kernels!.IncrementalOrder, () =>
lock (_gpu!.Value.KernelInitLock)
{
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>>>();
});
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>>>();
});
}
}
}
@@ -112,12 +129,15 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
return KernelContainer.InitOnce(ref _kernels!.Pow2ModOrder, () =>
lock (_gpu!.Value.KernelInitLock)
{
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>>>();
});
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>>>();
});
}
}
}
@@ -126,28 +146,68 @@ public sealed class GpuKernelLease
get
{
var accel = Accelerator;
return KernelContainer.InitOnce(ref _kernels!.SmallPrimeFactor, () =>
lock (_gpu!.Value.KernelInitLock)
{
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>>>();
});
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>>>();
});
}
}
}
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()
{
_stream?.Dispose();
_stream = null;
_gpu?.Dispose();
_gpu = null;
_releaseLimiter?.Invoke();
_releaseLimiter = null;
_kernels = null;
Pool.Enqueue(this);
Dispose(disposing: false);
}
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);
}
}
}
}
+3 -1
View File
@@ -280,7 +280,7 @@ public class GpuKernelPool
public static GpuKernelLease GetKernel(bool useGpuOrder)
{
var limiter = GpuPrimeWorkLimiter.Acquire();
var gpu = Rent();
var gpu = RentPreferred(preferCpu: !useGpuOrder);
var accelerator = gpu.Accelerator;
var kernels = GetKernels(accelerator);
return GpuKernelLease.Rent(limiter, gpu, kernels);
@@ -293,10 +293,12 @@ 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();
}
+31 -15
View File
@@ -3,42 +3,58 @@ 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 Lease Acquire()
public static IDisposable Acquire()
{
// TODO: Inline the pooled limiter guard from the GpuLimiterThroughputBenchmarks so prime-sieve
// GPU jobs stop allocating Lease instances on every acquisition.
// GPU jobs stop allocating Releaser instances on every acquisition.
var sem = _semaphore;
// sem.Wait();
return new Lease(sem);
sem.Wait();
return new Releaser(sem);
}
public static void SetLimit(int value)
{
if (value == _currentLimit)
{
return;
}
if (value < 1)
{
value = 1;
}
// 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;
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;
}
}
public sealed class Lease
private sealed class Releaser : IDisposable
{
private readonly SemaphoreSlim _sem;
private bool _disposed;
internal Lease(SemaphoreSlim sem)
public Releaser(SemaphoreSlim sem)
{
_sem = sem;
}
public void Dispose()
{
// _sem.Release();
if (_disposed)
{
return;
}
_sem.Release();
_disposed = true;
}
}
}
+13 -6
View File
@@ -6,14 +6,14 @@ public static class GpuWorkLimiter
private static readonly object _lock = new();
private static int _currentLimit = 1;
public static Lease Acquire()
public static IDisposable Acquire()
{
// TODO: Replace the per-call Lease allocation with the pooled struct-based guard from the
// TODO: Replace the per-call Releaser 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 Lease(sem);
return new Releaser(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 Lease instances.
// that may still be in use by existing Releaser 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,18 +40,25 @@ public static class GpuWorkLimiter
}
}
public sealed class Lease
private sealed class Releaser : IDisposable
{
private readonly SemaphoreSlim _sem;
private bool _disposed;
internal Lease(SemaphoreSlim sem)
public Releaser(SemaphoreSlim sem)
{
_sem = sem;
}
public void Dispose()
{
if (_disposed)
{
return;
}
_sem.Release();
_disposed = true;
}
}
}
@@ -1,97 +0,0 @@
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;
}
}
@@ -53,28 +53,15 @@ internal static class DivisorByDivisorKernels
int endIndex = offset + count;
if (hasCycle)
{
GpuUInt128 remainderValue = new(cycleRemainder);
for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++)
{
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;
remainderValue.AddMod(delta, cycleLength);
cycleRemainder = remainderValue.Low;
byte hit = 0;
if (cycleRemainder == 0UL)
@@ -1,73 +1,33 @@
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<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)
public static void SmallPrimeSieveKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<uint> smallPrimes, ArrayView<byte> results)
{
ulong n = numbers[index];
byte result = 1;
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;
}
// 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.
int length = (int)primes.Length;
long length = smallPrimes.Length;
for (int i = 0; i < length; i++)
{
ulong primeSquare = primeSquares[i];
ulong prime = smallPrimes[i];
ulong primeSquare = (ulong)prime * prime;
if (primeSquare > n)
{
break;
}
ulong prime = primes[i];
if (n % prime == 0UL)
{
result = 0;
@@ -78,32 +38,16 @@ internal static class PrimeTesterKernels
results[index] = result;
}
public static void HeuristicTrialDivisionKernel(
Index1D index,
ArrayView<int> resultFlag,
ulong n,
ulong maxDivisor,
HeuristicGpuDivisorTableKind tableKind,
HeuristicGpuDivisorTables tables,
HeuristicGpuCombinedDivisorTables combinedTables)
public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView<ulong> divisors, ulong n, ArrayView<byte> results)
{
byte nMod10 = (byte)(n % 10UL);
ArrayView1D<ulong, Stride1D.Dense> divisors = tableKind == HeuristicGpuDivisorTableKind.Combined
? combinedTables.SelectCombined(nMod10)
: tables.SelectDivisors(tableKind, nMod10);
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)
ulong divisor = divisors[index];
if (divisor <= 1UL)
{
resultFlag[0] = 1;
results[index] = 0;
return;
}
results[index] = n % divisor == 0UL ? (byte)1 : (byte)0;
}
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,6 +112,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
resources.HitIndexBuffer,
resources.Divisors,
resources.Exponents,
resources.FilteredDivisors,
resources.DivisorData,
resources.Offsets,
resources.Counts,
@@ -121,7 +122,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
out processedCount);
ReturnBatchResources(resources);
// Monitor.Exit(gpuLease.ExecutionLock);
Monitor.Exit(gpuLease.ExecutionLock);
ReturnAccelerator(gpuLease);
if (composite)
@@ -176,6 +177,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
MemoryBuffer1D<int, Stride1D.Dense> hitIndexBuffer,
ulong[] divisors,
ulong[] exponents,
ulong[] filteredDivisors,
GpuDivisorPartialData[] divisorData,
int[] offsets,
int[] counts,
@@ -217,23 +219,24 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
ulong remainingCount = maxK;
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);
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);
UInt128 currentDivisor128 = firstDivisor128;
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);
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));
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();
@@ -268,50 +271,37 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
chunkCount = (int)remainingCount;
}
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++)
{
bool passesSmallModuli = localRemainder3 != 0 && localRemainder7 != 0 && localRemainder11 != 0;
if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0)
{
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;
}
ulong candidate = (ulong)nextDivisor128;
divisorSpan[admissibleCount] = divisorValue;
divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue);
offsetSpan[admissibleCount] = admissibleCount;
countSpan[admissibleCount] = 1;
cycleSpan[admissibleCount] = divisorCycle;
exponentSpan[admissibleCount] = prime;
admissibleCount++;
if (DivisorGenerator.IsValidDivisor(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, decimalMask))
{
filteredStorage[filteredCount++] = candidate;
}
nextDivisor128 += twoP128;
localRemainder10 = AddMod10(localRemainder10, step10);
localRemainder10 = AddMod(localRemainder10, step10, 10);
localRemainder8 = AddMod8(localRemainder8, step8);
localRemainder3 = AddMod3(localRemainder3, step3);
localRemainder7 = AddMod7(localRemainder7, step7);
localRemainder11 = AddMod11(localRemainder11, step11);
localRemainder5 = AddMod(localRemainder5, step5, 5);
localRemainder3 = AddMod(localRemainder3, step3, 3);
localRemainder7 = AddMod(localRemainder7, step7, 7);
localRemainder11 = AddMod(localRemainder11, step11, 11);
}
remainder10 = localRemainder10;
remainder8 = localRemainder8;
remainder5 = localRemainder5;
remainder3 = localRemainder3;
remainder7 = localRemainder7;
remainder11 = localRemainder11;
@@ -323,6 +313,32 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
currentDivisor128 = nextDivisor128;
remainingCount -= (ulong)chunkCount;
if (filteredCount == 0)
{
continue;
}
int admissibleCount = 0;
for (int i = 0; i < filteredCount; i++)
{
ulong divisorValue = filteredStorage[i];
MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData);
if (divisorCycle != prime)
{
continue;
}
divisorSpan[admissibleCount] = divisorValue;
exponentSpan[admissibleCount] = prime;
divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue);
offsetSpan[admissibleCount] = admissibleCount;
countSpan[admissibleCount] = 1;
cycleSpan[admissibleCount] = divisorCycle;
admissibleCount++;
}
if (admissibleCount == 0)
{
continue;
@@ -366,62 +382,19 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static byte AddMod3(byte value, byte delta)
private static byte AddMod(byte value, byte delta, byte modulus)
{
const int Modulus = 3;
int sum = value + delta;
int result = value + delta;
if (sum >= Modulus)
if (result >= modulus)
{
sum -= Modulus;
result -= modulus;
}
return (byte)sum;
return (byte)result;
}
[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)
{
@@ -429,6 +402,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime)
{
// The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1,
@@ -452,25 +428,13 @@ 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 cappedLimit;
return divisorLimit;
}
ulong computedLimit = (1UL << (int)(prime - 1UL)) - 1UL;
if (prime < 10_000_000UL && computedLimit > cappedLimit)
{
computedLimit = cappedLimit;
}
return computedLimit < cappedLimit ? computedLimit : cappedLimit;
return computedLimit < divisorLimit ? computedLimit : divisorLimit;
}
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -488,17 +452,19 @@ 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.Rent();
_lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator;
_capacity = Math.Max(1, owner._gpuBatchSize);
}
internal void Reset()
{
_disposed = false;
}
private void EnsureExecutionResourcesLocked(int requiredCapacity)
@@ -545,6 +511,11 @@ 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)
@@ -580,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
}
// Monitor.Enter(_lease.ExecutionLock);
Monitor.Enter(_lease.ExecutionLock);
EnsureExecutionResourcesLocked(length);
@@ -613,7 +584,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
int sentinel = int.MaxValue;
firstHitView.CopyFromCPU(ref sentinel, 1);
var kernel = _kernel!;
var kernel = _kernel;
if (kernel is null)
{
Monitor.Exit(_lease.ExecutionLock);
throw new InvalidOperationException("GPU divisor kernel was not initialized.");
}
kernel(1, divisorView, offsetView, countView, exponentView, cycleView, hitView, firstHitView);
@@ -621,11 +597,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
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);
}
}
@@ -657,7 +639,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
return lease;
}
return GpuContextPool.Rent();
return GpuContextPool.RentPreferred(preferCpu: false);
}
private void ReturnAccelerator(GpuContextPool.GpuContextLease lease) => _acceleratorPool.Add(lease);
@@ -675,7 +657,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private void ReturnBatchResources(BatchResources resources) => _resourcePool.Add(resources);
private sealed class BatchResources
private sealed class BatchResources : IDisposable
{
private Accelerator? _accelerator;
private MemoryBuffer1D<GpuDivisorPartialData, Stride1D.Dense>? _divisorDataBuffer;
@@ -712,6 +694,8 @@ 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!;
@@ -786,6 +770,7 @@ 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);
@@ -824,12 +809,14 @@ 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.Rent();
var gpu = GpuContextPool.RentPreferred(preferCpu: false);
var accelerator = gpu.Accelerator;
var kernel = GetKernel(accelerator);
var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D<byte>(1));
@@ -17,6 +17,7 @@ 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
@@ -46,82 +47,88 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
int currentSize;
int i;
UInt128 q = UInt128.Zero;
while (kStart <= maxK && Volatile.Read(ref isPrime))
try
{
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)
while (kStart <= maxK && Volatile.Read(ref isPrime))
{
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++)
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)
{
if (orders[i] != 0UL)
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++)
{
Volatile.Write(ref isPrime, false);
break;
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)
else
{
if (orders[i] != 0UL && q.IsPrimeCandidate())
q = twoP * kStart + UInt128.One;
for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
{
Volatile.Write(ref isPrime, false);
if (orders[i] != 0UL && q.IsPrimeCandidate())
{
Volatile.Write(ref isPrime, false);
}
}
}
}
kStart += (UInt128)currentSize;
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,6 +38,16 @@ 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:
@@ -68,7 +78,7 @@ public class MersenneNumberLucasLehmerGpuTester
}
else
{
var gpu = GpuContextPool.Rent();
var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu);
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.
@@ -100,7 +110,7 @@ public class MersenneNumberLucasLehmerGpuTester
}
}
var gpu = GpuContextPool.Rent();
var gpu = GpuContextPool.RentPreferred(preferCpu: false);
var accelerator = gpu.Accelerator;
var kernel = GetBatchKernel(accelerator);
@@ -142,17 +152,23 @@ public class MersenneNumberLucasLehmerGpuTester
}
GpuUInt128[] buffer = ArrayPool<GpuUInt128>.Shared.Rent(count);
ComputeResidues(exponents, buffer.AsSpan(0, count));
for (int i = 0; i < count; i++)
try
{
results[i] = buffer[i].IsZero;
ComputeResidues(exponents, buffer.AsSpan(0, count));
for (int i = 0; i < count; i++)
{
results[i] = buffer[i].IsZero;
}
}
finally
{
ArrayPool<GpuUInt128>.Shared.Return(buffer);
}
ArrayPool<GpuUInt128>.Shared.Return(buffer);
}
private bool IsMersennePrimeNtt(ulong exponent, GpuUInt128 nttMod, GpuUInt128 primitiveRoot, bool runOnGpu)
{
var gpu = GpuContextPool.Rent();
var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu);
var accelerator = gpu.Accelerator;
var addKernel = GetAddSmallKernel(accelerator);
var subKernel = GetSubSmallKernel(accelerator);
@@ -4,18 +4,19 @@ 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 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 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
var kernel = _kernelType switch
{
@@ -25,41 +26,48 @@ 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);
UInt128 remaining;
int currentSize;
int found = 0;
while (kStart <= maxK && Volatile.Read(ref isPrime))
try
{
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);
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);
// 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;
}
// 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);
kStart += (UInt128)currentSize;
}
stream.Synchronize();
foundBuffer.View.CopyToCPU(ref found, 1);
if (found != 0)
{
Volatile.Write(ref isPrime, false);
break;
}
foundBuffer.Dispose();
gpuLease.Dispose();
}
kStart += (UInt128)currentSize;
}
}
finally
{
foundBuffer.Dispose();
execution.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.Rent();
return GpuContextPool.RentPreferred(preferCpu: !_useGpuOrder);
}
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
private sealed class ResidueResources : IDisposable
{
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
private sealed class SquareCacheEntry : IDisposable
{
public MemoryBuffer1D<GpuUInt128, Stride1D.Dense> Buffer { get; }
public GpuUInt128 Root { get; }
@@ -194,20 +194,21 @@ 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();
@@ -241,12 +242,13 @@ 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;
}
@@ -360,11 +362,12 @@ 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
@@ -464,6 +467,7 @@ internal static partial class PrimeOrderGpuHeuristics
stackProductBuffer.Dispose();
resultBuffer.Dispose();
statusBuffer.Dispose();
execution.Dispose();
}
}
@@ -522,6 +526,7 @@ 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);
@@ -546,6 +551,7 @@ 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;
}
@@ -553,6 +559,7 @@ internal static partial class PrimeOrderGpuHeuristics
remainderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), exponents.Length);
exponentBuffer.Dispose();
remainderBuffer.Dispose();
execution.Dispose();
lease.Dispose();
return true;
}
@@ -580,6 +587,7 @@ internal static partial class PrimeOrderGpuHeuristics
try
{
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var kernel = GetPow2ModWideKernel(accelerator);
@@ -620,6 +628,7 @@ internal static partial class PrimeOrderGpuHeuristics
exponentBuffer.Dispose();
remainderBuffer.Dispose();
execution.Dispose();
return true;
}
catch (Exception)
@@ -1,717 +0,0 @@
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);
}
}
}
+19 -20
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 ulong[] GroupADivisorsStorage;
private static readonly int[] GroupADivisorsStorage;
private static readonly ulong[] GroupADivisorSquaresStorage;
private static readonly ulong[] GroupBDivisorsEnding1Storage;
private static readonly int[] GroupBDivisorsEnding1Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding1Storage;
private static readonly ulong[] GroupBDivisorsEnding3Storage;
private static readonly int[] GroupBDivisorsEnding3Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding3Storage;
private static readonly ulong[] GroupBDivisorsEnding7Storage;
private static readonly int[] GroupBDivisorsEnding7Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding7Storage;
private static readonly ulong[] GroupBDivisorsEnding9Storage;
private static readonly int[] GroupBDivisorsEnding9Storage;
private static readonly ulong[] GroupBDivisorSquaresEnding9Storage;
static HeuristicPrimeSieves()
@@ -46,46 +46,45 @@ public static class HeuristicPrimeSieves
// Intentionally left blank. Accessing this method forces the static constructor to run.
}
public static ReadOnlySpan<ulong> GroupADivisors => GroupADivisorsStorage;
public static ReadOnlySpan<int> GroupADivisors => GroupADivisorsStorage;
public static ReadOnlySpan<ulong> GroupADivisorSquares => GroupADivisorSquaresStorage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding1 => GroupBDivisorsEnding1Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding1 => GroupBDivisorsEnding1Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding1 => GroupBDivisorSquaresEnding1Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding3 => GroupBDivisorsEnding3Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding3 => GroupBDivisorsEnding3Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding3 => GroupBDivisorSquaresEnding3Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding7 => GroupBDivisorsEnding7Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding7 => GroupBDivisorsEnding7Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding7 => GroupBDivisorSquaresEnding7Storage;
public static ReadOnlySpan<ulong> GroupBDivisorsEnding9 => GroupBDivisorsEnding9Storage;
public static ReadOnlySpan<int> GroupBDivisorsEnding9 => GroupBDivisorsEnding9Storage;
public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding9 => GroupBDivisorSquaresEnding9Storage;
private static (ulong[] divisors, ulong[] squares) BuildGroupADivisors()
private static (int[] divisors, ulong[] squares) BuildGroupADivisors()
{
var divisors = new ulong[SieveLength];
var divisors = new int[SieveLength];
var squares = new ulong[SieveLength];
int index = 0;
for (; index < GroupAConstants.Length && index < divisors.Length; index++)
{
ulong value = (ulong)GroupAConstants[index];
int value = GroupAConstants[index];
divisors[index] = value;
squares[index] = value * value;
squares[index] = (ulong)value * (ulong)value;
}
int candidate = 23;
int incrementIndex = 0;
while (index < divisors.Length)
{
ulong value = (ulong)candidate;
divisors[index] = value;
squares[index] = value * value;
divisors[index] = candidate;
squares[index] = (ulong)candidate * (ulong)candidate;
candidate += GroupAIncrementPattern[incrementIndex];
incrementIndex ^= 1;
@@ -95,9 +94,9 @@ public static class HeuristicPrimeSieves
return (divisors, squares);
}
private static (ulong[] divisors, ulong[] squares) BuildGroupBResidueSieve(ReadOnlySpan<ushort> residues)
private static (int[] divisors, ulong[] squares) BuildGroupBResidueSieve(ReadOnlySpan<ushort> residues)
{
var divisors = new ulong[SieveLength];
var divisors = new int[SieveLength];
var squares = new ulong[SieveLength];
int index = 0;
@@ -113,7 +112,7 @@ public static class HeuristicPrimeSieves
continue;
}
divisors[index] = candidate;
divisors[index] = (int)candidate;
squares[index] = candidate * candidate;
index++;
}
File diff suppressed because it is too large. Load diff
@@ -16,10 +16,8 @@ public interface IMersenneNumberDivisorByDivisorTester
void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues);
public interface IDivisorScanSession
public interface IDivisorScanSession : IDisposable
{
void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan<ulong> primes, Span<byte> hits);
void Dispose();
}
}
+76 -82
View File
@@ -4,7 +4,6 @@ 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;
@@ -137,31 +136,20 @@ 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)
{
int read = 0;
Span<byte> span = buffer;
while (read < span.Length)
ulong d, c;
try
{
int chunk = stream.Read(span[read..]);
if (chunk == 0)
{
if (read == 0)
{
yield break;
}
throw new EndOfStreamException("Incomplete divisor-cycle entry.");
}
read += chunk;
d = reader.ReadUInt64();
c = reader.ReadUInt64();
}
catch (EndOfStreamException)
{
yield break;
}
ulong divisor = BinaryPrimitives.ReadUInt64LittleEndian(span[..8]);
ulong cycle = BinaryPrimitives.ReadUInt64LittleEndian(span[8..]);
yield return (divisor, cycle);
yield return (d, c);
}
}
@@ -391,8 +379,7 @@ public class MersenneDivisorCycles
for (int i = 0; i < smallLength; i++)
{
ulong prime = smallPrimes[i];
ulong primeSquared = smallPrimesSquared[i];
if (primeSquared > remaining)
if (smallPrimesSquared[i] > remaining)
{
break;
}
@@ -411,9 +398,7 @@ public class MersenneDivisorCycles
// Atomic.Add(ref _tryFactorIntoCountsInternalHits, 1UL);
// Console.WriteLine($"MersenneDivisorCycles.TryFactorIntoCountsInternal hits {Volatile.Read(ref _tryFactorIntoCountsInternalHits)}");
// HeuristicPrimeTester primeTester = _primeTester ??= new();
if (PrimeTester.IsPrime(remaining))
// if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None))
if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None))
{
AddFactor(counts, remaining);
return true;
@@ -453,14 +438,26 @@ public class MersenneDivisorCycles
private static void AddFactor(Dictionary<ulong, int> counts, ulong factor)
{
ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists);
entry = exists ? entry + 1 : 1;
if (counts.TryGetValue(factor, out int existing))
{
counts[factor] = existing + 1;
}
else
{
counts[factor] = 1;
}
}
private static void AddFactor(Dictionary<ulong, int> counts, ulong factor, int multiplicity)
{
ref int entry = ref CollectionsMarshal.GetValueRefOrAddDefault(counts, factor, out bool exists);
entry = exists ? entry + multiplicity : multiplicity;
if (counts.TryGetValue(factor, out int existing))
{
counts[factor] = existing + multiplicity;
}
else
{
counts[factor] = multiplicity;
}
}
private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
@@ -538,8 +535,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 writeIndex = multiplicity;
while (writeIndex > 0)
int actual = 0;
while (actual < multiplicity)
{
if (working < prime)
{
@@ -552,20 +549,26 @@ public class MersenneDivisorCycles
break;
}
writeIndex--;
candidateBuffer[writeIndex] = reduced;
candidateBuffer[actual] = reduced;
working = reduced;
actual++;
}
int actual = multiplicity - writeIndex;
if (actual == 0)
{
return;
}
Span<ulong> candidates = candidateBuffer.Slice(writeIndex, actual);
Span<ulong> candidates = candidateBuffer.Slice(0, 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]);
@@ -595,28 +598,24 @@ public class MersenneDivisorCycles
// return 2UL;
// }
ulong modulus = n;
ulong nMinus1 = modulus - 1UL;
ulong nMinus2 = modulus - 2UL;
while (true)
{
ulong c = (NextRandomUInt64() % nMinus1) + 1UL;
ulong x = (NextRandomUInt64() % nMinus2) + 2UL;
ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL;
ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL;
ulong y = x;
ulong d = 1UL;
while (d == 1UL)
{
x = AdvancePolynomial(x, c, modulus);
y = AdvancePolynomial(y, c, modulus);
y = AdvancePolynomial(y, c, modulus);
x = AdvancePolynomial(x, c, n);
y = AdvancePolynomial(y, c, n);
y = AdvancePolynomial(y, c, n);
ulong diff = x > y ? x - y : y - x;
d = BinaryGcd(diff, modulus);
d = BinaryGcd(diff, n);
}
if (d != modulus)
if (d != n)
{
return d;
}
@@ -917,22 +916,19 @@ public class MersenneDivisorCycles
return CalculateCycleLengthFallback(divisor);
}
// [ThreadStatic]
// private static HeuristicPrimeTester? _primeTester;
internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false)
{
if ((divisor & (divisor - 1UL)) == 0UL)
{
cycleLength = 1UL;
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 (divisor <= 3UL)
{
cycleLength = CalculateCycleLengthFallback(divisor);
return true;
}
// if (!skipPrimeOrderHeuristic)
// {
@@ -940,27 +936,25 @@ public class MersenneDivisorCycles
// Console.WriteLine($"MersenneDivisorCycles.CalculateCycleLengthFallback hits {Volatile.Read(ref _tryCalculateCycleLengthHeuristicHits)}");
// }
// 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;
}
}
// 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;
}
}
cycleLength = 0UL;
return false;
}
cycleLength = 0UL;
return false;
}
private static ulong CalculateCycleLengthFallback(ulong divisor)
{
@@ -1,6 +1,5 @@
using System.Buffers;
using System.Threading;
using System.Threading.Tasks;
using System.Collections.Concurrent;
using System.Globalization;
using System.Runtime.InteropServices;
@@ -168,10 +167,14 @@ public static class MersenneNumberDivisorByDivisorTester
if (primesToTest.Count == 0)
{
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);
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.");
}
return;
}
@@ -242,20 +245,19 @@ public static class MersenneNumberDivisorByDivisorTester
void ProcessPrime(ulong prime)
{
Console.WriteLine($"Processing {prime}");
Console.WriteLine($"Task started {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($"Finished processing {prime}");
Console.WriteLine($"Task finished {prime}");
}
if (workerCount == 1)
@@ -267,7 +269,12 @@ public static class MersenneNumberDivisorByDivisorTester
}
else
{
TaskScheduler scheduler = UnboundedTaskScheduler.Instance;
var scheduler = UnboundedTaskScheduler.Instance;
ParallelOptions options = new()
{
MaxDegreeOfParallelism = workerCount,
TaskScheduler = scheduler
};
int totalCount = filteredPrimes.Count;
int partitionSize = chunkSize < 1 ? 1 : chunkSize;
@@ -276,41 +283,13 @@ public static class MersenneNumberDivisorByDivisorTester
partitionSize = totalCount;
}
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)
Parallel.ForEach(Partitioner.Create(0, totalCount, partitionSize), options, range =>
{
int rangeStart = start;
int rangeEnd = rangeStart + partitionSize;
if (rangeEnd > totalCount)
for (int index = range.Item1; index < range.Item2; index++)
{
rangeEnd = totalCount;
ProcessPrime(filteredPrimes[index]);
}
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();
});
}
}
}
+44 -37
View File
@@ -147,54 +147,61 @@ 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!;
// 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)
try
{
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++)
// 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)
{
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)
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++)
{
order = qs[offset + i].CalculateOrder();
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;
}
}
if (useOrderCache)
{
OrderCache[qs[offset + i]] = order;
}
offset += count;
orderBuffer.Dispose();
qBuffer.Dispose();
}
offset += count;
orderBuffer.Dispose();
qBuffer.Dispose();
}
gpuLease.Dispose();
ArrayPool<UInt128>.Shared.Return(qs);
finally
{
execution.Dispose();
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)
if (q < 2UL || (q & 1UL) == 0)
{
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.IsPrime(q))
if (!PrimeTester.IsPrimeCpu(q, ct))
{
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.IsPrime(m))
if (PrimeTester.IsPrimeCpu(m, ct))
{
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.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
+222 -150
View File
@@ -299,16 +299,22 @@ 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();
int last = actual - 1;
evaluations[last] = stepper.InitializeCpuIsUnity(candidates[last]);
for (int j = last - 1; j >= 0; j--)
evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]);
for (int j = 1; j < actual; j++)
{
evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]);
}
for (int j = 0; j < actual; j++)
for (int j = actual - 1; j >= 0; j--)
{
if (!evaluations[j])
{
@@ -458,16 +464,16 @@ internal static partial class PrimeOrderCalculator
}
Span<ulong> candidates = buffer[..actual];
candidates.Reverse();
stepper.Reset();
int last = actual - 1;
if (stepper.InitializeCpuIsUnity(candidates[last]))
if (stepper.InitializeCpuIsUnity(candidates[0]))
{
return true;
}
for (int j = last - 1; j >= 0; j--)
for (int j = 1; j < actual; j++)
{
if (stepper.ComputeNextIsUnity(candidates[j]))
{
@@ -688,16 +694,13 @@ 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, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3);
CandidateKey keyY = BuildKey(y, previous, previousGroup, hasPrevious, threshold1, threshold2, threshold3);
CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious);
CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious);
int compare = keyX.Primary.CompareTo(keyY.Primary);
if (compare != 0)
{
@@ -714,9 +717,9 @@ internal static partial class PrimeOrderCalculator
});
}
private static CandidateKey BuildKey(ulong value, ulong previous, int previousGroup, bool hasPrevious, ulong threshold1, ulong threshold2, ulong threshold3)
private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious)
{
int group = GetGroup(value, threshold1, threshold2, threshold3);
int group = GetGroup(value, prime);
if (group == 0)
{
return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue);
@@ -766,18 +769,21 @@ internal static partial class PrimeOrderCalculator
return groupOffset + (isGe ? 0 : 1);
}
private static int GetGroup(ulong value, ulong threshold1, ulong threshold2, ulong threshold3)
private static int GetGroup(ulong value, ulong prime)
{
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;
@@ -822,11 +828,9 @@ internal static partial class PrimeOrderCalculator
}
FactorEntry factor = factors[index];
int factorExponent = factor.Exponent;
ulong primeFactor = factor.Value;
ulong contribution = 1UL;
ulong contributionLimit = factorExponent == 0 ? 0UL : order / primeFactor;
for (int exponent = 0; exponent <= factorExponent; exponent++)
for (int exponent = 0; exponent <= factor.Exponent; exponent++)
{
ulong nextDivisor = divisorProduct * contribution;
if (nextDivisor > order)
@@ -840,12 +844,12 @@ internal static partial class PrimeOrderCalculator
return;
}
if (exponent == factorExponent)
if (exponent == factor.Exponent)
{
break;
}
if (contribution > contributionLimit)
if (contribution > order / primeFactor)
{
break;
}
@@ -969,45 +973,24 @@ internal static partial class PrimeOrderCalculator
}
Span<ulong> candidates = buffer[..actual];
candidates.Reverse();
stepper.Reset();
int last = actual - 1;
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)
if (powUsed >= powBudget && powBudget > 0)
{
return true;
}
powUsed++;
if (stepper.InitializeCpuIsUnity(candidates[last]))
if (stepper.InitializeCpuIsUnity(candidates[0]))
{
return true;
}
for (int j = last - 1; j >= 0; j--)
for (int j = 1; j < actual; j++)
{
if (powUsed >= powBudget)
if (powUsed >= powBudget && powBudget > 0)
{
return true;
}
@@ -1057,13 +1040,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;
@@ -1078,20 +1061,45 @@ 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)
{
PopulateSmallPrimeFactorsCpu(
value,
limit,
ref primeSlots,
ref exponentSlots,
ref primeSlotArray,
ref exponentSlotArray,
out factorCount,
out remaining);
// 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;
}
}
@@ -1125,13 +1133,12 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorHits, 1UL);
// Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}");
// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
bool isPrime = PrimeTester.IsPrime(composite);
bool isPrime = PrimeTester.IsPrimeGpu(composite);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
if (isPrime)
{
AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1);
AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
continue;
}
@@ -1193,8 +1200,7 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorPendingHits, 1UL);
// Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}");
// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
bool isPrime = PrimeTester.IsPrime(composite);
bool isPrime = PrimeTester.IsPrimeGpu(composite);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
entry = entry.WithPrimality(isPrime);
@@ -1203,7 +1209,7 @@ internal static partial class PrimeOrderCalculator
if (entry.IsPrime)
{
AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1);
AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
}
// composite is never smaller on the execution path
// else if (composite > 1UL)
@@ -1231,22 +1237,39 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
cofactorIsPrime = PrimeTester.IsPrime(cofactor);
// cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
cofactorIsPrime = PrimeTester.IsPrimeCpu(cofactor, CancellationToken.None);
// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
}
ArrayPool<FactorEntry> pool = ThreadStaticPools.FactorEntryPool;
bool temp;
if (factorCount == 0)
if (useDictionary)
{
if (cofactor == value)
temp = counts is null || counts.Count == 0;
if (temp && cofactor == value)
{
result = PartialFactorResult.Rent(null, cofactor, false, 0, cofactorIsPrime);
goto ReturnResult;
}
result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0, cofactorIsPrime);
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);
goto ReturnResult;
}
@@ -1286,10 +1309,14 @@ internal static partial class PrimeOrderCalculator
result = PartialFactorResult.Rent(array, cofactor, temp, arrayIndex, cofactorIsPrime);
ReturnResult:
if (primeSlotArray is not null)
// ThreadStaticPools.UlongPool.Return(primeSlotsArray);
// ThreadStaticPools.IntPool.Return(exponentSlotsArray);
if (counts is not null)
{
ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false);
ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false);
counts.Clear();
ThreadStaticPools.ReturnUlongIntDictionary(counts);
}
if (compositeStack is not null)
@@ -1303,28 +1330,35 @@ internal static partial class PrimeOrderCalculator
return result;
}
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)
private static bool TryPopulateSmallPrimeFactorsCpu(
ulong value,
uint limit,
Span<ulong> primeTargets,
Span<int> exponentTargets,
out int factorCount,
out ulong remaining)
{
factorCount = 0;
ulong remainingLocal = value;
remaining = value;
if (value <= 1UL)
{
remaining = remainingLocal;
return;
}
// 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;
// }
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++)
@@ -1343,13 +1377,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;
@@ -1361,13 +1395,10 @@ internal static partial class PrimeOrderCalculator
continue;
}
EnsureCollectorCapacity(
ref primeTargets,
ref exponentTargets,
ref primeArray,
ref exponentArray,
factorCount,
factorCount + 1);
if (factorCount >= capacity)
{
throw new InvalidOperationException($"Capacity is smaller than factor count");
}
primeTargets[factorCount] = primeValue;
exponentTargets[factorCount] = exponent;
@@ -1375,6 +1406,7 @@ internal static partial class PrimeOrderCalculator
}
remaining = remainingLocal;
return true;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -1442,6 +1474,46 @@ 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)
@@ -1564,18 +1636,11 @@ internal static partial class PrimeOrderCalculator
}
}
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)
private static bool TryAppendFactor(Span<ulong> primes, Span<int> exponents, ref int count, ulong prime, int exponent)
{
if (exponent <= 0)
{
return;
return true;
}
for (int i = 0; i < count; i++)
@@ -1583,60 +1648,70 @@ internal static partial class PrimeOrderCalculator
if (primes[i] == prime)
{
exponents[i] += exponent;
return;
return true;
}
}
EnsureCollectorCapacity(
ref primes,
ref exponents,
ref primeArray,
ref exponentArray,
count,
count + 1);
if (count >= primes.Length)
{
return false;
}
primes[count] = prime;
exponents[count] = exponent;
count++;
return true;
}
private static void EnsureCollectorCapacity(
ref Span<ulong> primes,
ref Span<int> exponents,
ref ulong[]? primeArray,
ref int[]? exponentArray,
int count,
int requiredCapacity)
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)
{
if (requiredCapacity <= primes.Length)
if (exponent <= 0)
{
return;
}
int newCapacity = primes.Length << 1;
if (newCapacity < requiredCapacity)
if (useDictionary)
{
newCapacity = requiredCapacity;
AddFactor(counts!, prime, exponent);
return;
}
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)
if (TryAppendFactor(primes, exponents, ref count, prime, exponent))
{
ThreadStaticPools.UlongPool.Return(primeArray, clearArray: false);
ThreadStaticPools.IntPool.Return(exponentArray!, clearArray: false);
return;
}
primeArray = newPrimeArray;
exponentArray = newExponentArray;
primes = newPrimeSpan;
exponents = newExponentSpan;
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;
}
}
private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData)
@@ -1683,11 +1758,9 @@ 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)
{
@@ -1698,11 +1771,11 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _factorCompletelyHits, 1UL);
// Console.WriteLine($"Factor completely hits {Volatile.Read(ref _factorCompletelyHits)}");
// HeuristicPrimeTester tester = _tester ??= new();
// PrimeTester primeTester = PrimeTester.Exclusive;
// bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
bool isPrime = PrimeTester.IsPrime(value);
// bool isPrime = PrimeTester.IsPrimeGpu(value);
// bool isPrime = PrimeTester.IsPrimeCpu(value, CancellationToken.None);
bool isPrime = PrimeTester.IsPrimeCpu(value, CancellationToken.None);
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
if (!knownComposite && isPrime)
{
@@ -1717,8 +1790,8 @@ internal static partial class PrimeOrderCalculator
// Console.WriteLine($"Factor completely after PollardRho hits {Volatile.Read(ref _factorCompletelyPollardRhoHits)}");
// isPrime = primeTester.HeuristicIsPrimeGpu(factor);
isPrime = PrimeTester.IsPrime(factor);
// isPrime = PrimeTester.IsPrimeGpu(factor);
// isPrime = PrimeTester.IsPrimeCpu(factor, CancellationToken.None);
isPrime = PrimeTester.IsPrimeCpu(factor, CancellationToken.None);
// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
if (isPrime)
{
@@ -1754,8 +1827,7 @@ internal static partial class PrimeOrderCalculator
// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
isPrime = PrimeTester.IsPrime(quotient);
// isPrime = PrimeTester.IsPrimeGpu(quotient);
isPrime = PrimeTester.IsPrimeCpu(quotient, CancellationToken.None);
// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
}
+71 -48
View File
@@ -9,63 +9,86 @@ 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);
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);
// primeBuffer.Clear();
// exponentBuffer.Clear();
remaining = value;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
try
{
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
try
{
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator);
SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots);
SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator);
SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots);
// scratch.Clear();
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];
counts.Add(primeValue, exponent);
}
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;
// }
lease.Dispose();
ThreadStaticPools.UlongPool.Return(primeBufferArray);
ThreadStaticPools.IntPool.Return(exponentBufferArray);
return true;
}
counts.Add(primeValue, exponent);
}
return true;
}
finally
{
execution.Dispose();
lease.Dispose();
}
}
finally
{
ThreadStaticPools.UlongPool.Return(primeBufferArray);
ThreadStaticPools.IntPool.Return(exponentBufferArray);
}
}
}
@@ -185,49 +185,55 @@ 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);
Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
factorSpan.CopyTo(buffer);
if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor))
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
{
buffer[length] = new FactorEntry128(factors.Cofactor, 1);
length++;
}
buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value));
Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
factorSpan.CopyTo(buffer);
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 (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor))
{
if ((order % primeFactor) == UInt128.Zero)
{
UInt128 reduced = order / primeFactor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
order = reduced;
continue;
}
}
break;
buffer[length] = new FactorEntry128(factors.Cofactor, 1);
length++;
}
}
pool.Return(tempArray, clearArray: false);
return order;
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)
{
UInt128 reduced = order / primeFactor;
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
order = reduced;
continue;
}
}
break;
}
}
return order;
}
finally
{
pool.Return(tempArray, clearArray: false);
}
}
private static bool TryConfirmOrderWide(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
@@ -950,7 +956,7 @@ internal static partial class PrimeOrderCalculator
if (value <= ulong.MaxValue)
{
// return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None);
return PrimeTester.IsPrime((ulong)value);
return PrimeTester.IsPrimeCpu((ulong)value, CancellationToken.None);
}
if ((value & UInt128.One) == UInt128.Zero)
File diff suppressed because it is too large. Load diff
+146 -6
View File
@@ -1,11 +1,151 @@
using System;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
public static class PrimesGenerator
{
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;
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,
};
}
}
+34 -150
View File
@@ -6,7 +6,7 @@ using System.Threading.Tasks;
namespace PerfectNumbers.Core;
internal sealed class TaskThreadPool
internal sealed class TaskThreadPool : IDisposable
{
private readonly struct WorkItem
{
@@ -19,15 +19,10 @@ internal sealed class TaskThreadPool
}
private readonly ConcurrentQueue<WorkItem> _pendingTaskQueue = new();
private readonly HashSet<Thread> _runningThreads = new();
private readonly object _lock = new();
private readonly BlockingCollection<WorkItem> _pendingTasks;
private readonly Thread[] _threads;
private readonly Action<Task> _taskExecutor;
private readonly int _maximumConcurrency;
private readonly SemaphoreSlim _workAvailable = new(0);
private readonly CancellationTokenSource _disposeCancellation = new();
private int _pendingTaskCount;
private int _currentThreads;
private int _threadId;
private int _disposed;
public TaskThreadPool(int minimumThreads, Action<Task> taskExecutor)
{
@@ -37,28 +32,23 @@ internal sealed class TaskThreadPool
}
_taskExecutor = taskExecutor ?? throw new ArgumentNullException(nameof(taskExecutor));
_maximumConcurrency = minimumThreads;
_pendingTasks = new BlockingCollection<WorkItem>(_pendingTaskQueue);
_threads = new Thread[minimumThreads];
List<Thread> threadsToStart = new(minimumThreads);
lock (_lock)
for (int i = 0; i < minimumThreads; i++)
{
for (int i = 0; i < minimumThreads; i++)
Thread thread = new(WorkerLoop)
{
Thread thread = CreateWorkerThread();
_runningThreads.Add(thread);
_currentThreads++;
threadsToStart.Add(thread);
}
}
IsBackground = true,
Name = $"{nameof(TaskThreadPool)}-{i + 1}"
};
foreach (Thread thread in threadsToStart)
{
_threads[i] = thread;
thread.Start();
}
}
public int ThreadCount => _maximumConcurrency;
public int ThreadCount => _threads.Length;
public void Queue(Task task)
{
@@ -67,29 +57,22 @@ internal sealed class TaskThreadPool
throw new ArgumentNullException(nameof(task));
}
List<Thread>? threadsToStart;
lock (_lock)
if (Volatile.Read(ref _disposed) != 0)
{
if (_disposeCancellation.IsCancellationRequested)
{
throw new ObjectDisposedException(nameof(TaskThreadPool));
}
_pendingTaskQueue.Enqueue(new WorkItem(task));
_pendingTaskCount++;
threadsToStart = EnsureWorkerCapacityLocked(null);
throw new ObjectDisposedException(nameof(TaskThreadPool));
}
_workAvailable.Release();
if (threadsToStart is not null)
try
{
foreach (Thread thread in threadsToStart)
{
thread.Start();
}
_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);
}
}
@@ -114,132 +97,33 @@ internal sealed class TaskThreadPool
public void Dispose()
{
Thread[] threadsToJoin;
lock (_lock)
if (Interlocked.Exchange(ref _disposed, 1) != 0)
{
_disposeCancellation.Cancel();
threadsToJoin = new Thread[_runningThreads.Count];
_runningThreads.CopyTo(threadsToJoin);
return;
}
if (threadsToJoin.Length > 0)
{
_workAvailable.Release(threadsToJoin.Length);
}
_pendingTasks.CompleteAdding();
foreach (Thread thread in threadsToJoin)
foreach (Thread thread in _threads)
{
thread.Join();
}
_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;
}
}
_pendingTasks.Dispose();
}
private void WorkerLoop()
{
try
{
if (!TryDequeueWork(out WorkItem workItem))
foreach (WorkItem workItem in _pendingTasks.GetConsumingEnumerable())
{
return;
}
_taskExecutor(workItem.Task);
}
finally
{
CompleteWorker();
}
}
private void CompleteWorker()
{
List<Thread>? threadsToStart = null;
lock (_lock)
{
_runningThreads.Remove(Thread.CurrentThread);
_currentThreads--;
if (!_disposeCancellation.IsCancellationRequested)
{
threadsToStart = EnsureWorkerCapacityLocked(null);
_taskExecutor(workItem.Task);
}
}
if (threadsToStart is not null)
catch (ObjectDisposedException)
{
foreach (Thread thread in threadsToStart)
{
thread.Start();
}
// The collection was disposed while the worker was waiting. Exit gracefully.
}
}
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
public sealed class TextFileWriter : IDisposable
{
private readonly FileStream _stream;
private StreamWriter? _textWriter;
+37
View File
@@ -260,6 +260,43 @@ 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.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(result, modulus, nPrime);
result = result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(multiplier, modulus, nPrime);
result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
index = windowStart - 1;
}
@@ -75,7 +75,7 @@ public static partial class ULongExtensions
if (!keepMontgomery)
{
result = result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
result = result.MontgomeryMultiply(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.MontgomeryMultiplyCpu(square, modulus, nPrime);
oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
}
}
+9 -47
View File
@@ -46,45 +46,6 @@ 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)
{
@@ -243,15 +204,14 @@ public static partial class ULongExtensions
ulong modulus = divisor.Modulus;
ulong nPrime = divisor.NPrime;
ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent);
return MontgomeryMultiplyGpu(montgomeryResult, 1UL, modulus, nPrime);
return MontgomeryMultiply(montgomeryResult, 1UL, modulus, nPrime);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorData divisor, ulong exponent)
{
ulong modulus = divisor.Modulus;
// Reuse the GPU-compatible MontgomeryMultiply helper highlighted in MontgomeryMultiplyBenchmarks so the accelerator path
// mirrors the fastest CPU reduction without relying on UInt128 intrinsics.
// Reuse the optimized MontgomeryMultiply extension highlighted in MontgomeryMultiplyBenchmarks so the GPU path matches the fastest CPU reduction.
// We shouldn't hit it in production code. If we do, we're doing something wrong.
// if (exponent == 0UL)
// {
@@ -273,7 +233,7 @@ public static partial class ULongExtensions
while (index >= 0)
{
ulong currentBit = (exponent >> index) & 1UL;
ulong squared = MontgomeryMultiplyGpu(result, result, modulus, nPrime);
ulong squared = MontgomeryMultiply(result, result, modulus, nPrime);
bool processWindow = currentBit != 0UL;
result = processWindow ? result : squared;
@@ -293,11 +253,11 @@ public static partial class ULongExtensions
int windowLength = processWindow ? index - windowStart + 1 : 0;
for (int square = 0; square < windowLength; square++)
{
result = MontgomeryMultiplyGpu(result, result, modulus, nPrime);
result = MontgomeryMultiply(result, result, modulus, nPrime);
}
result = processWindow
? MontgomeryMultiplyGpu(
? MontgomeryMultiply(
result,
ComputeMontgomeryOddPowerGpu(
(exponent >> windowStart) & ((1UL << windowLength) - 1UL),
@@ -370,7 +330,7 @@ public static partial class ULongExtensions
while (remaining != 0UL)
{
ulong bit = remaining & 1UL;
ulong multiplied = MontgomeryMultiplyGpu(power, baseValue, modulus, nPrime);
ulong multiplied = MontgomeryMultiply(power, baseValue, modulus, nPrime);
ulong mask = (ulong)-(long)bit;
power = (power & ~mask) | (multiplied & mask);
@@ -380,7 +340,7 @@ public static partial class ULongExtensions
break;
}
baseValue = MontgomeryMultiplyGpu(baseValue, baseValue, modulus, nPrime);
baseValue = MontgomeryMultiply(baseValue, baseValue, modulus, nPrime);
}
return power;
@@ -457,6 +417,7 @@ 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.
@@ -483,6 +444,7 @@ 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 MontgomeryMultiplyCpu(this ulong a, ulong b, ulong modulus, ulong nPrime)
public static ulong MontgomeryMultiply(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.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
}
remainingExponent >>= 1;
@@ -590,10 +590,10 @@ public static partial class ULongExtensions
break;
}
baseVal = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
}
return keepMontgomery ? result : result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime);
}
+1 -1
View File
@@ -131,7 +131,7 @@ 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 `maxDivisorSquare`/`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 `sqrtLimit`/`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
View File
@@ -1,8 +0,0 @@
# 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.