MontgomeryDivisorDataPool.Rent / .Return removed as we no longer pool Montgomery divisor data. UInt128Extensions.MontgomeryDivisorDataPool variables optimized.

This commit is contained in:
mlemanczyk committed 2025-12-16 02:27:29 +01:00
1 parent b7773adad8
commit 3c80221105
27 files changed
+257 -473

No files matched your search

@@ -500,7 +500,7 @@ public class Pow2ModBenchmarks
Console.WriteLine($"Calculating small cycle {smallModulus}"); Console.WriteLine($"Calculating small cycle {smallModulus}");
smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
smallExponents[i] = NextSmallExponent(random); smallExponents[i] = NextSmallExponent(random);
smallCycles[i] = MersenneDivisorCycles.CalculateCycleLengthCpu(smallModulus, MontgomeryDivisorDataPool.Shared.FromModulus(smallModulus)); smallCycles[i] = MersenneDivisorCycles.CalculateCycleLengthCpu(smallModulus, MontgomeryDivisorData.FromModulus(smallModulus));
(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(random); (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(random);
Console.WriteLine($"Calculating large cycle {largeModulus}"); Console.WriteLine($"Calculating large cycle {largeModulus}");
@@ -613,7 +613,7 @@ public class Pow2ModBenchmarks
private static ulong CalculateCycleLengthWithHeuristicsCpu(ulong modulus, ref ulong? previousPrimeOrder) private static ulong CalculateCycleLengthWithHeuristicsCpu(ulong modulus, ref ulong? previousPrimeOrder)
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
return MersenneDivisorCycles.CalculateCycleLengthCpu(modulus, divisorData); return MersenneDivisorCycles.CalculateCycleLengthCpu(modulus, divisorData);
@@ -680,7 +680,7 @@ public class Pow2ModBenchmarks
private ulong CalculateCycleLengthWithHeuristicsCpu(ulong modulus, ulong fallbackOrder = 0UL) private ulong CalculateCycleLengthWithHeuristicsCpu(ulong modulus, ulong fallbackOrder = 0UL)
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLengthCpu(modulus, divisorData); return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLengthCpu(modulus, divisorData);
@@ -709,7 +709,7 @@ public class Pow2ModBenchmarks
private ulong CalculateCycleLengthWithHeuristicsGpu(PrimeOrderCalculatorAccelerator gpu, ulong modulus, ulong fallbackOrder = 0UL) private ulong CalculateCycleLengthWithHeuristicsGpu(PrimeOrderCalculatorAccelerator gpu, ulong modulus, ulong fallbackOrder = 0UL)
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, modulus, divisorData); return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, modulus, divisorData);
@@ -39,7 +39,7 @@ public class Pow2MontgomeryModBenchmarks
ulong smallModulus = NextSmallOddModulus(); ulong smallModulus = NextSmallOddModulus();
_smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
_smallExponents[i] = NextSmallExponent(); _smallExponents[i] = NextSmallExponent();
_smallCycles[i] = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, smallModulus, MontgomeryDivisorDataPool.Shared.FromModulus(smallModulus)); _smallCycles[i] = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, smallModulus, MontgomeryDivisorData.FromModulus(smallModulus));
Console.WriteLine($"Calculating large cycle {i + 1} with heuristics"); Console.WriteLine($"Calculating large cycle {i + 1} with heuristics");
(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(); (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
@@ -66,7 +66,7 @@ internal static class Pow2MontgomeryModWindowedBenchmarkInputsProvider
{ {
Pow2MontgomeryModWindowedBenchmarkSeed seed = seeds[i]; Pow2MontgomeryModWindowedBenchmarkSeed seed = seeds[i];
ulong modulus = 2UL * seed.Multiplier * seed.Exponent + 1UL; ulong modulus = 2UL * seed.Multiplier * seed.Exponent + 1UL;
MontgomeryDivisorData divisor = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus);
ulong cycleLength = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, modulus, divisor); ulong cycleLength = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, modulus, divisor);
ulong reducedExponent = seed.Exponent % cycleLength; ulong reducedExponent = seed.Exponent % cycleLength;
cases[i] = new Pow2MontgomeryModWindowedBenchmarkCase( cases[i] = new Pow2MontgomeryModWindowedBenchmarkCase(
@@ -5,7 +5,7 @@ using PerfectNumbers.Core.Gpu;
using PerfectNumbers.Core.Gpu.Accelerators; using PerfectNumbers.Core.Gpu.Accelerators;
using Xunit; using Xunit;
namespace PerfectNumbers.Core.Tests; namespace PerfectNumbers.Core.Tests.Gpu;
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
public class PrimeOrderGpuHeuristicsTests public class PrimeOrderGpuHeuristicsTests
@@ -25,7 +25,7 @@ public class PrimeOrderGpuHeuristicsTests
try try
{ {
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, 1UL, prime, out ulong remainder, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, 1UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow); status.Should().Be(GpuPow2ModStatus.Overflow);
@@ -43,7 +43,7 @@ public class PrimeOrderGpuHeuristicsTests
{ {
const ulong exponent = 138_000_001UL; const ulong exponent = 138_000_001UL;
const ulong prime = 138_000_001UL; const ulong prime = 138_000_001UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, prime, out ulong remainder, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, prime, out ulong remainder, divisorData);
@@ -62,7 +62,7 @@ public class PrimeOrderGpuHeuristicsTests
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
try try
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, 1UL, prime, out ulong remainder, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, 1UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow); status.Should().Be(GpuPow2ModStatus.Overflow);
@@ -85,7 +85,7 @@ public class PrimeOrderGpuHeuristicsTests
try try
{ {
const ulong prime = 193UL; const ulong prime = 193UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, 16UL, prime, out ulong remainder, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, 16UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow); status.Should().Be(GpuPow2ModStatus.Overflow);
@@ -113,7 +113,7 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] remainders = new ulong[2]; ulong[] remainders = new ulong[2];
ulong[] exponents = { 1UL, 2UL }; ulong[] exponents = { 1UL, 2UL };
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow); status.Should().Be(GpuPow2ModStatus.Overflow);
@@ -133,7 +133,7 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] remainders = { 123UL, 456UL }; ulong[] remainders = { 123UL, 456UL };
ulong[] exponents = { 138_000_001UL, 138_000_011UL }; ulong[] exponents = { 138_000_001UL, 138_000_011UL };
const ulong prime = 138_000_001UL; const ulong prime = 138_000_001UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData);
@@ -162,7 +162,7 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] exponents = { 1UL, 2UL }; ulong[] exponents = { 1UL, 2UL };
ulong[] remainders = { 5UL, 6UL }; ulong[] remainders = { 5UL, 6UL };
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow); status.Should().Be(GpuPow2ModStatus.Overflow);
@@ -189,7 +189,7 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] remainders = { 9UL, 11UL }; ulong[] remainders = { 9UL, 11UL };
const ulong prime = 193UL; const ulong prime = 193UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow); status.Should().Be(GpuPow2ModStatus.Overflow);
@@ -397,7 +397,7 @@ public class PrimeOrderGpuHeuristicsTests
{ {
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
ulong result = PrimeOrderCalculator.CalculateGpu(gpu, prime, null, divisorData, PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); ulong result = PrimeOrderCalculator.CalculateGpu(gpu, prime, null, divisorData, PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault);
result.Should().Be(3UL); result.Should().Be(3UL);
@@ -18,7 +18,7 @@ public class ULongExtensionsGpuTests
[InlineData(65535UL, 2147483647UL)] [InlineData(65535UL, 2147483647UL)]
public void Pow2ModWindowedGpu_matches_montgomery_cpu_for_known_values(ulong exponent, ulong modulus) public void Pow2ModWindowedGpu_matches_montgomery_cpu_for_known_values(ulong exponent, ulong modulus)
{ {
MontgomeryDivisorData divisor = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus);
ulong expected = exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisor); ulong expected = exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisor);
exponent.Pow2ModWindowedGpu(modulus).Should().Be(expected); exponent.Pow2ModWindowedGpu(modulus).Should().Be(expected);
@@ -52,7 +52,7 @@ public class ULongExtensionsGpuTests
foreach (ulong modulus in moduli) foreach (ulong modulus in moduli)
{ {
MontgomeryDivisorData divisor = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus);
foreach (ulong exponent in exponents) foreach (ulong exponent in exponents)
{ {
ulong expected = exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisor); ulong expected = exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisor);
@@ -40,7 +40,7 @@ public class MersenneDivisorCyclesGenerateGpuTests
foreach (var pair in pairs) foreach (var pair in pairs)
{ {
expected.Should().Contain(pair.divisor); expected.Should().Contain(pair.divisor);
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(pair.divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor);
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLengthCpu(pair.divisor, divisorData)); pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLengthCpu(pair.divisor, divisorData));
} }
} }
@@ -90,7 +90,7 @@ public class MersenneDivisorCyclesGenerateGpuTests
foreach (var pair in pairs) foreach (var pair in pairs)
{ {
uniqueDivisors.Add(pair.divisor).Should().BeTrue(); uniqueDivisors.Add(pair.divisor).Should().BeTrue();
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(pair.divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor);
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLengthCpu(pair.divisor, divisorData)); pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLengthCpu(pair.divisor, divisorData));
} }
@@ -22,7 +22,7 @@ public class MersenneDivisorCyclesTests
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
try try
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
MersenneDivisorCycles.CalculateCycleLengthCpu(divisor, divisorData).Should().Be(expected); MersenneDivisorCycles.CalculateCycleLengthCpu(divisor, divisorData).Should().Be(expected);
MersenneDivisorCycles.CalculateCycleLengthGpu(divisor).Should().Be(expected); MersenneDivisorCycles.CalculateCycleLengthGpu(divisor).Should().Be(expected);
} }
@@ -49,17 +49,17 @@ public class MersenneNumberDivisorCpuTesterTests
{ {
ulong[] primes = [5UL, 7UL, 11UL, 13UL]; ulong[] primes = [5UL, 7UL, 11UL, 13UL];
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLengthCpu(23UL, MontgomeryDivisorDataPool.Shared.FromModulus(23UL)); ulong cycle23 = MersenneDivisorCycles.CalculateCycleLengthCpu(23UL, MontgomeryDivisorData.FromModulus(23UL));
session.CheckDivisor(23UL, MontgomeryDivisorDataPool.Shared.FromModulus(23UL), cycle23, primes).Should().BeTrue(); session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes).Should().BeTrue();
ulong[] primesWithoutHit = [5UL, 7UL, 13UL]; ulong[] primesWithoutHit = [5UL, 7UL, 13UL];
session.CheckDivisor(23UL, MontgomeryDivisorDataPool.Shared.FromModulus(23UL), cycle23, primesWithoutHit).Should().BeFalse(); session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primesWithoutHit).Should().BeFalse();
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLengthCpu(31UL, MontgomeryDivisorDataPool.Shared.FromModulus(31UL)); ulong cycle31 = MersenneDivisorCycles.CalculateCycleLengthCpu(31UL, MontgomeryDivisorData.FromModulus(31UL));
session.CheckDivisor(31UL, MontgomeryDivisorDataPool.Shared.FromModulus(31UL), cycle31, primes).Should().BeTrue(); session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes).Should().BeTrue();
ulong[] primesBeforeHit = [7UL, 11UL, 13UL]; ulong[] primesBeforeHit = [7UL, 11UL, 13UL];
session.CheckDivisor(31UL, MontgomeryDivisorDataPool.Shared.FromModulus(31UL), cycle31, primesBeforeHit).Should().BeFalse(); session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primesBeforeHit).Should().BeFalse();
} }
finally finally
{ {
@@ -81,7 +81,7 @@ public class MersenneNumberDivisorCpuTesterTests
{ {
ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL]; ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL];
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(7UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(7UL, divisorData); ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents).Should().BeTrue(); session.CheckDivisor(7UL, divisorData, cycle, exponents).Should().BeTrue();
@@ -108,7 +108,7 @@ public class MersenneNumberDivisorCpuTesterTests
{ {
ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL]; ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL];
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(11UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(11UL, divisorData); ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents).Should().BeTrue(); session.CheckDivisor(11UL, divisorData, cycle, exponents).Should().BeTrue();
@@ -231,10 +231,10 @@ public class MersenneNumberDivisorGpuTesterTests
{ {
ulong[] primes = [31UL, 37UL, 41UL]; ulong[] primes = [31UL, 37UL, 41UL];
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), 0UL, primes).Should().BeTrue(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes).Should().BeTrue();
ulong[] primesWithoutHit = [31UL, 41UL]; ulong[] primesWithoutHit = [31UL, 41UL];
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), 0UL, primesWithoutHit).Should().BeFalse(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primesWithoutHit).Should().BeFalse();
} }
finally finally
{ {
@@ -256,17 +256,17 @@ public class MersenneNumberDivisorGpuTesterTests
{ {
ulong[] primes = [31UL, 37UL, 41UL, 43UL]; ulong[] primes = [31UL, 37UL, 41UL, 43UL];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL)); ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primes).Should().BeTrue(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue();
ulong[] primesBeforeHit = [31UL]; ulong[] primesBeforeHit = [31UL];
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse();
ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 13367UL, MontgomeryDivisorDataPool.Shared.FromModulus(13367UL)); ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 13367UL, MontgomeryDivisorData.FromModulus(13367UL));
session.CheckDivisor(13367UL, MontgomeryDivisorDataPool.Shared.FromModulus(13367UL), cycle13367, primes).Should().BeTrue(); session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes).Should().BeTrue();
ulong[] primesWithoutMatchingExponent = [31UL, 37UL, 43UL]; ulong[] primesWithoutMatchingExponent = [31UL, 37UL, 43UL];
session.CheckDivisor(13367UL, MontgomeryDivisorDataPool.Shared.FromModulus(13367UL), cycle13367, primesWithoutMatchingExponent).Should().BeFalse(); session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primesWithoutMatchingExponent).Should().BeFalse();
} }
finally finally
{ {
@@ -292,11 +292,11 @@ public class MersenneNumberDivisorGpuTesterTests
{ {
ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL]; ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL)); ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primes).Should().BeTrue(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue();
ulong[] primesBeforeHit = [31UL]; ulong[] primesBeforeHit = [31UL];
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse();
} }
finally finally
{ {
@@ -321,7 +321,7 @@ public class MersenneNumberDivisorGpuTesterTests
{ {
ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL]; ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL];
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(7UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 7UL, divisorData); ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents).Should().BeTrue(); session.CheckDivisor(7UL, divisorData, cycle, exponents).Should().BeTrue();
@@ -352,7 +352,7 @@ public class MersenneNumberDivisorGpuTesterTests
{ {
ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL]; ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL];
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(11UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 11UL, divisorData); ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents).Should().BeTrue(); session.CheckDivisor(11UL, divisorData, cycle, exponents).Should().BeTrue();
@@ -384,7 +384,7 @@ public class MersenneNumberDivisorGpuTesterTests
var patchedTable = new List<(ulong Divisor, ulong Cycle)>(); var patchedTable = new List<(ulong Divisor, ulong Cycle)>();
ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1];
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(191UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(191UL);
patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 191UL, divisorData); patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 191UL, divisorData);
tableField.SetValue(cycles, patchedTable); tableField.SetValue(cycles, patchedTable);
@@ -428,11 +428,11 @@ public class MersenneNumberDivisorGpuTesterTests
{ {
ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL]; ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL)); ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primes).Should().BeTrue(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue();
ulong[] primesBeforeHit = [31UL]; ulong[] primesBeforeHit = [31UL];
session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse();
} }
finally finally
{ {
@@ -49,7 +49,7 @@ public class PrimeOrderCalculatorTests
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
MontgomeryDivisorDataPool.Shared.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
@@ -69,14 +69,14 @@ public class PrimeOrderCalculatorTests
gpu, gpu,
prime, prime,
previousOrder, previousOrder,
MontgomeryDivisorDataPool.Shared.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault);
ulong strict = PrimeOrderCalculator.CalculateGpu( ulong strict = PrimeOrderCalculator.CalculateGpu(
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
MontgomeryDivisorDataPool.Shared.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.StrictDefault); PrimeOrderCalculator.PrimeOrderCalculatorConfig.StrictDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
@@ -95,7 +95,7 @@ public class PrimeOrderCalculatorTests
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
MontgomeryDivisorDataPool.Shared.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
@@ -117,7 +117,7 @@ public class PrimeOrderCalculatorTests
gpu, gpu,
prime: 13UL, prime: 13UL,
previousOrder: null, previousOrder: null,
MontgomeryDivisorDataPool.Shared.FromModulus(13UL), MontgomeryDivisorData.FromModulus(13UL),
config); config);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
@@ -137,7 +137,7 @@ public class PrimeOrderCalculatorTests
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
try try
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(239UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
ulong heuristic = PrimeOrderCalculator.CalculateGpu( ulong heuristic = PrimeOrderCalculator.CalculateGpu(
gpu, gpu,
prime: 239UL, prime: 239UL,
@@ -5,20 +5,20 @@ namespace PerfectNumbers.Core.Cpu;
internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{ {
private Func<ulong, ulong> _getCycleLength = orderDevice switch private Func<ulong, MontgomeryDivisorData, ulong> _getCycleLength = orderDevice switch
{ {
ComputationDevice.Gpu => divisor => DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor), ComputationDevice.Gpu => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor, divisorData),
ComputationDevice.Hybrid => divisor => DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor), ComputationDevice.Hybrid => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor, divisorData),
_ => static (divisor) => DivisorCycleCache.Shared.GetCycleLengthCpu(divisor), _ => static (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthCpu(divisor, divisorData),
}; };
public void Configure(PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) public void Configure(PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice)
{ {
_getCycleLength = orderDevice switch _getCycleLength = orderDevice switch
{ {
ComputationDevice.Gpu => divisor => DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor), ComputationDevice.Gpu => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor, divisorData),
ComputationDevice.Hybrid => divisor => DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor), ComputationDevice.Hybrid => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor, divisorData),
_ => static (divisor) => DivisorCycleCache.Shared.GetCycleLengthCpu(divisor), _ => static (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthCpu(divisor, divisorData),
}; };
} }
@@ -45,7 +45,7 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera
if (divisorCycle == 0UL) if (divisorCycle == 0UL)
{ {
divisorCycle = _getCycleLength(divisor); divisorCycle = _getCycleLength(divisor, divisorData);
// This never kicks off in production code // This never kicks off in production code
// if (divisorCycle == 0UL) // if (divisorCycle == 0UL)
// { // {
@@ -6,18 +6,29 @@ using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core.Cpu; namespace PerfectNumbers.Core.Cpu;
public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester public sealed class MersenneNumberDivisorByDivisorCpuTester(ComputationDevice orderDevice = ComputationDevice.Hybrid) : IMersenneNumberDivisorByDivisorTester
{ {
private BigInteger _divisorLimit; private BigInteger _divisorLimit;
private int _batchSize = 1_024; private int _batchSize = 1_024;
private BigInteger _minK = BigInteger.One; private BigInteger _minK = BigInteger.One;
private readonly ComputationDevice _orderDevice; private readonly ComputationDevice OrderDevice = orderDevice;
private readonly TryCycleLengthDelegate _tryCalculateCycleLengthForExponent; private readonly TryCycleLengthDelegate _tryCalculateCycleLengthForExponent = orderDevice switch
private readonly CalculateCycleLengthDelegate _calculateCycleLength; {
ComputationDevice.Gpu => TryCalculateCycleLengthForExponentGpu,
ComputationDevice.Hybrid => TryCalculateCycleLengthForExponentHybrid,
_ => TryCalculateCycleLengthForExponentCpu,
};
private readonly CalculateCycleLengthDelegate _calculateCycleLength = orderDevice switch
{
ComputationDevice.Gpu => CalculateCycleLengthGpu,
ComputationDevice.Hybrid => CalculateCycleLengthHybrid,
_ => CalculateCycleLengthCpu,
};
private string? _stateFilePath; private string? _stateFilePath;
private int _stateCounter; private int _stateCounter;
private BigInteger _lastSavedK; private BigInteger _lastSavedK;
public ComputationDevice OrderDevice => _orderDevice;
private delegate bool TryCycleLengthDelegate( private delegate bool TryCycleLengthDelegate(
ulong divisor, ulong divisor,
@@ -43,24 +54,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
public GpuUInt128 Limit; public GpuUInt128 Limit;
} }
public MersenneNumberDivisorByDivisorCpuTester(ComputationDevice orderDevice = ComputationDevice.Hybrid)
{
_orderDevice = orderDevice;
_tryCalculateCycleLengthForExponent = orderDevice switch
{
ComputationDevice.Gpu => TryCalculateCycleLengthForExponentGpu,
ComputationDevice.Hybrid => TryCalculateCycleLengthForExponentHybrid,
_ => TryCalculateCycleLengthForExponentCpu,
};
_calculateCycleLength = orderDevice switch
{
ComputationDevice.Gpu => CalculateCycleLengthGpu,
ComputationDevice.Hybrid => CalculateCycleLengthHybrid,
_ => CalculateCycleLengthCpu,
};
}
public int BatchSize public int BatchSize
{ {
get => _batchSize; get => _batchSize;
@@ -157,7 +150,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
} }
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) => ThreadStaticPools.RentMersenneCpuDivisorSession(gpu, _orderDevice); public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) => ThreadStaticPools.RentMersenneCpuDivisorSession(gpu, OrderDevice);
private bool CheckDivisors( private bool CheckDivisors(
PrimeOrderCalculatorAccelerator gpu, PrimeOrderCalculatorAccelerator gpu,
@@ -340,13 +333,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
ulong divisorHigh = divisor.High; ulong divisorHigh = divisor.High;
var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor);
var divisorPool = MontgomeryDivisorDataPool.Shared;
while (divisor.CompareTo(limit) <= 0) while (divisor.CompareTo(limit) <= 0)
{ {
if (residueStepper.IsAdmissible()) if (residueStepper.IsAdmissible())
{ {
ulong candidate = divisor.Low; ulong candidate = divisor.Low;
MontgomeryDivisorData divisorData = divisorPool.FromModulus(candidate); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
ulong divisorCycle; 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. // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs.
if (!_tryCalculateCycleLengthForExponent( if (!_tryCalculateCycleLengthForExponent(
@@ -370,7 +362,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
divisorCycle = computedCycle; divisorCycle = computedCycle;
} }
divisorPool.Return(divisorData);
RecordState(currentK); RecordState(currentK);
if (divisorCycle == prime) if (divisorCycle == prime)
{ {
@@ -553,12 +544,11 @@ private bool CheckDivisors64(
out ulong foundDivisor) out ulong foundDivisor)
{ {
bool canAdvance = step <= limit; bool canAdvance = step <= limit;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
if (!canAdvance) if (!canAdvance)
{ {
if (divisor <= limit && residueStepper.IsAdmissible()) if (divisor <= limit && residueStepper.IsAdmissible())
{ {
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
ulong divisorCycle; ulong divisorCycle;
if (!_tryCalculateCycleLengthForExponent( if (!_tryCalculateCycleLengthForExponent(
divisor, divisor,
@@ -579,8 +569,6 @@ private bool CheckDivisors64(
divisorCycle = computedCycle; divisorCycle = computedCycle;
} }
divisorPool.Return(divisorData);
if (divisorCycle == prime) if (divisorCycle == prime)
{ {
foundDivisor = divisor; foundDivisor = divisor;
@@ -611,7 +599,7 @@ private bool CheckDivisors64(
{ {
if (residueStepper.IsAdmissible()) if (residueStepper.IsAdmissible())
{ {
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
ulong divisorCycle; ulong divisorCycle;
if (!_tryCalculateCycleLengthForExponent( if (!_tryCalculateCycleLengthForExponent(
divisor, divisor,
@@ -632,7 +620,6 @@ private bool CheckDivisors64(
divisorCycle = computedCycle; divisorCycle = computedCycle;
} }
divisorPool.Return(divisorData);
RecordState(currentK); RecordState(currentK);
if (divisorCycle == prime) if (divisorCycle == prime)
{ {
+62 -149
View File
@@ -19,22 +19,18 @@ public sealed class DivisorCycleCache
private static int _sharedDivisorCyclesBatchSize = GpuConstants.GpuCycleStepsPerInvocation; private static int _sharedDivisorCyclesBatchSize = GpuConstants.GpuCycleStepsPerInvocation;
private const int CycleCacheTrackingLimit = 500_000; private const int CycleCacheTrackingLimit = 500_000;
private static readonly HashSet<ulong> CycleCacheTrackedDivisors = new(); private static readonly HashSet<ulong> CycleCacheTrackedDivisors = [];
private static readonly object CycleCacheTrackingLock = new(); private static readonly object CycleCacheTrackingLock = new();
private static ulong CycleCacheHitCount;
private static bool CycleCacheTrackingDisabled;
public static void SetDivisorCyclesBatchSize(int divisorCyclesBatchSize) public static void SetDivisorCyclesBatchSize(int divisorCyclesBatchSize)
{ {
_sharedDivisorCyclesBatchSize = Math.Max(1, divisorCyclesBatchSize); _sharedDivisorCyclesBatchSize = Math.Max(1, divisorCyclesBatchSize);
} }
public static DivisorCycleCache Shared { get; } = new DivisorCycleCache(_sharedDivisorCyclesBatchSize); public static readonly DivisorCycleCache Shared = new(_sharedDivisorCyclesBatchSize);
public int PreferredBatchSize => _divisorCyclesBatchSize; public int PreferredBatchSize => _divisorCyclesBatchSize;
// private static readonly Accelerator[] _accelerators;
private DivisorCycleCache(int divisorCyclesBatchSize) private DivisorCycleCache(int divisorCyclesBatchSize)
{ {
_divisorCyclesBatchSize = Math.Max(1, divisorCyclesBatchSize); _divisorCyclesBatchSize = Math.Max(1, divisorCyclesBatchSize);
@@ -51,54 +47,28 @@ public sealed class DivisorCycleCache
_snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); _snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
} }
public ulong GetCycleLengthCpu(ulong divisor, bool skipPrimeOrderHeuristic = false) [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public ulong GetCycleLengthCpu(ulong divisor, in MontgomeryDivisorData divisorData)
{ {
if (!skipPrimeOrderHeuristic)
{
// TODO: Implement array-free version of GetCycleLengthsCpu so that we don't pay the array costs for single checks.
Span<ulong> result = stackalloc ulong[1];
ReadOnlySpan<ulong> singleDivisor = stackalloc ulong[1] { divisor };
GetCycleLengthsCpu(singleDivisor, result);
return result[0];
}
// The skipPrimeOrderHeuristic path is only used for divisors produced by the CPU by-divisor scanner,
// which always generates values greater than one. Keep this guard disabled here to avoid repeating the
// validation already performed by the general GetCycleLengths path.
// if (divisor <= 1UL)
// {
// throw new InvalidDataException("Divisor must be > 1");
// }
ulong[] snapshot = _snapshot; ulong[] snapshot = _snapshot;
ulong cycleLength;
if (divisor < (ulong)snapshot.Length) if (divisor < (ulong)snapshot.Length)
{ {
cycleLength = snapshot[divisor]; ulong cached = snapshot[divisor];
if (cycleLength == 0UL) ArgumentOutOfRangeException.ThrowIfZero(cached);
{
throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
}
return cycleLength; return cached;
} }
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; return MersenneDivisorCycles.CalculateCycleLengthCpu(divisor, divisorData);
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
cycleLength = MersenneDivisorCycles.CalculateCycleLengthCpu(divisor, divisorData, skipPrimeOrderHeuristic: true);
divisorPool.Return(divisorData);
return cycleLength;
} }
public ulong GetCycleLengthHybrid(PrimeOrderCalculatorAccelerator gpu, ulong divisor, bool skipPrimeOrderHeuristic = false) [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public ulong GetCycleLengthHybrid(PrimeOrderCalculatorAccelerator gpu, ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false)
{ {
if (!skipPrimeOrderHeuristic) if (!skipPrimeOrderHeuristic)
{ {
// TODO: Implement array-free version of GetCycleLengthsHybrid so that we don't pay the array costs for single checks. return GetCycleLengthHybrid(gpu, divisor, divisorData);
Span<ulong> result = stackalloc ulong[1];
ReadOnlySpan<ulong> singleDivisor = stackalloc ulong[1] { divisor };
GetCycleLengthsHybrid(gpu, singleDivisor, result);
return result[0];
} }
// The skipPrimeOrderHeuristic path is only used for divisors produced by the CPU by-divisor scanner, // The skipPrimeOrderHeuristic path is only used for divisors produced by the CPU by-divisor scanner,
@@ -114,22 +84,17 @@ public sealed class DivisorCycleCache
if (divisor < (ulong)snapshot.Length) if (divisor < (ulong)snapshot.Length)
{ {
cycleLength = snapshot[divisor]; cycleLength = snapshot[divisor];
if (cycleLength == 0UL) ArgumentOutOfRangeException.ThrowIfZero(cycleLength);
{
throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
}
return cycleLength; return cycleLength;
} }
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
cycleLength = MersenneDivisorCycles.CalculateCycleLengthHybrid(gpu, divisor, divisorData, skipPrimeOrderHeuristic: true); cycleLength = MersenneDivisorCycles.CalculateCycleLengthHybrid(gpu, divisor, divisorData, skipPrimeOrderHeuristic: true);
divisorPool.Return(divisorData);
return cycleLength; return cycleLength;
} }
public ulong GetCycleLengthGpu(PrimeOrderCalculatorAccelerator gpu, ulong divisor, bool skipPrimeOrderHeuristic = false) [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public ulong GetCycleLengthGpu(PrimeOrderCalculatorAccelerator gpu, ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false)
{ {
if (!skipPrimeOrderHeuristic) if (!skipPrimeOrderHeuristic)
{ {
@@ -152,18 +117,12 @@ public sealed class DivisorCycleCache
if (divisor < (ulong)snapshot.Length) if (divisor < (ulong)snapshot.Length)
{ {
cycleLength = snapshot[divisor]; cycleLength = snapshot[divisor];
if (cycleLength == 0UL) ArgumentOutOfRangeException.ThrowIfZero(cycleLength);
{
throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
}
return cycleLength; return cycleLength;
} }
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor);
cycleLength = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, divisor, divisorData, skipPrimeOrderHeuristic: true); cycleLength = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, divisor, divisorData, skipPrimeOrderHeuristic: true);
divisorPool.Return(divisorData);
return cycleLength; return cycleLength;
} }
@@ -181,18 +140,11 @@ public sealed class DivisorCycleCache
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
ulong divisor = divisors[i]; ulong divisor = divisors[i];
if (divisor <= 1UL)
{
throw new InvalidDataException("Divisor must be > 1");
}
if (divisor < (ulong)snapshot.Length) if (divisor < (ulong)snapshot.Length)
{ {
ulong cached = snapshot[divisor]; ulong cached = snapshot[divisor];
if (cached == 0UL) ArgumentOutOfRangeException.ThrowIfZero(cached);
{
throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
}
cycles[i] = cached; cycles[i] = cached;
} }
@@ -214,6 +166,22 @@ public sealed class DivisorCycleCache
} }
} }
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public ulong GetCycleLengthHybrid(PrimeOrderCalculatorAccelerator gpu, ulong divisor, in MontgomeryDivisorData divisorData)
{
ulong[] snapshot = _snapshot;
if (divisor < (ulong)snapshot.Length)
{
ulong cached = snapshot[divisor];
ArgumentOutOfRangeException.ThrowIfZero(cached);
return cached;
}
return MersenneDivisorCycles.CalculateCycleLengthHybrid(gpu, divisor, divisorData);
}
public void GetCycleLengthsHybrid(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> divisors, Span<ulong> cycles) public void GetCycleLengthsHybrid(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> divisors, Span<ulong> cycles)
{ {
ulong[] snapshot = _snapshot; ulong[] snapshot = _snapshot;
@@ -228,18 +196,11 @@ public sealed class DivisorCycleCache
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
ulong divisor = divisors[i]; ulong divisor = divisors[i];
if (divisor <= 1UL)
{
throw new InvalidDataException("Divisor must be > 1");
}
if (divisor < (ulong)snapshot.Length) if (divisor < (ulong)snapshot.Length)
{ {
ulong cached = snapshot[divisor]; ulong cached = snapshot[divisor];
if (cached == 0UL) ArgumentOutOfRangeException.ThrowIfZero(cached);
{
throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
}
cycles[i] = cached; cycles[i] = cached;
} }
@@ -275,18 +236,11 @@ public sealed class DivisorCycleCache
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
ulong divisor = divisors[i]; ulong divisor = divisors[i];
if (divisor <= 1UL)
{
throw new InvalidDataException("Divisor must be > 1");
}
if (divisor < (ulong)snapshot.Length) if (divisor < (ulong)snapshot.Length)
{ {
ulong cached = snapshot[divisor]; ulong cached = snapshot[divisor];
if (cached == 0UL) ArgumentOutOfRangeException.ThrowIfZero(cached);
{
throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
}
cycles[i] = cached; cycles[i] = cached;
} }
@@ -308,68 +262,27 @@ public sealed class DivisorCycleCache
} }
} }
private static void ObserveCycleCacheDivisor(ulong divisor)
{
bool logHit = false;
ulong hitNumber = 0UL;
bool trackingDisabledNow = false;
lock (CycleCacheTrackingLock)
{
if (CycleCacheTrackingDisabled)
{
return;
}
if (!CycleCacheTrackedDivisors.Add(divisor))
{
CycleCacheHitCount++;
hitNumber = CycleCacheHitCount;
logHit = true;
}
else if (CycleCacheTrackedDivisors.Count >= CycleCacheTrackingLimit)
{
CycleCacheTrackingDisabled = true;
CycleCacheTrackedDivisors.Clear();
trackingDisabledNow = true;
}
}
if (logHit)
{
Console.WriteLine($"Cycle cache hit for divisor {divisor} ({hitNumber})");
}
else if (trackingDisabledNow)
{
Console.WriteLine("Cycle cache hit tracking disabled after exceeding the tracking limit.");
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ComputeCyclesCpu(ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices) private static void ComputeCyclesCpu(ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices)
{ {
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
for (int i = 0; i < indices.Length; i++) for (int i = 0; i < indices.Length; i++)
{ {
int targetIndex = indices[i]; int targetIndex = indices[i];
ulong divisor = divisors[targetIndex]; ulong divisor = divisors[targetIndex];
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLengthCpu(divisor, divisorData); cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLengthCpu(divisor, divisorData);
divisorPool.Return(divisorData);
} }
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ComputeCyclesHybrid(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices) private static void ComputeCyclesHybrid(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices)
{ {
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
for (int i = 0; i < indices.Length; i++) for (int i = 0; i < indices.Length; i++)
{ {
int targetIndex = indices[i]; int targetIndex = indices[i];
ulong divisor = divisors[targetIndex]; ulong divisor = divisors[targetIndex];
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLengthHybrid(gpu, divisor, divisorData); cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLengthHybrid(gpu, divisor, divisorData);
divisorPool.Return(divisorData);
} }
} }
@@ -408,44 +321,44 @@ public sealed class DivisorCycleCache
private const int DefaultCapacity = 128; private const int DefaultCapacity = 128;
[ThreadStatic] [ThreadStatic]
private static MemoryBuffer1D<ulong, Stride1D.Dense>? divisorBuffer; private static MemoryBuffer1D<ulong, Stride1D.Dense>? _divisorBuffer;
[ThreadStatic] [ThreadStatic]
private static MemoryBuffer1D<ulong, Stride1D.Dense>? powBuffer; private static MemoryBuffer1D<ulong, Stride1D.Dense>? _powBuffer;
[ThreadStatic] [ThreadStatic]
private static MemoryBuffer1D<ulong, Stride1D.Dense>? orderBuffer; private static MemoryBuffer1D<ulong, Stride1D.Dense>? _orderBuffer;
[ThreadStatic] [ThreadStatic]
private static MemoryBuffer1D<ulong, Stride1D.Dense>? resultBuffer; private static MemoryBuffer1D<ulong, Stride1D.Dense>? _resultBuffer;
[ThreadStatic] [ThreadStatic]
private static MemoryBuffer1D<byte, Stride1D.Dense>? statusBuffer; private static MemoryBuffer1D<byte, Stride1D.Dense>? _statusBuffer;
private void EnsureCapacity(int requiredCapacity) private void EnsureCapacity(int requiredCapacity)
{ {
if (requiredCapacity <= (divisorBuffer?.Length ?? 0)) if (_divisorBuffer is not null && requiredCapacity <= _divisorBuffer.Length)
{ {
return; return;
} }
if (divisorBuffer != null) if (_divisorBuffer != null)
{ {
divisorBuffer.Dispose(); _divisorBuffer.Dispose();
powBuffer!.Dispose(); _powBuffer!.Dispose();
orderBuffer!.Dispose(); _orderBuffer!.Dispose();
resultBuffer!.Dispose(); _resultBuffer!.Dispose();
statusBuffer!.Dispose(); _statusBuffer!.Dispose();
} }
Accelerator accelerator = _accelerator; Accelerator accelerator = _accelerator;
// lock (accelerator) // lock (accelerator)
{ {
divisorBuffer = accelerator.Allocate1D<ulong>(requiredCapacity); _divisorBuffer = accelerator.Allocate1D<ulong>(requiredCapacity);
powBuffer = accelerator.Allocate1D<ulong>(requiredCapacity); _powBuffer = accelerator.Allocate1D<ulong>(requiredCapacity);
orderBuffer = accelerator.Allocate1D<ulong>(requiredCapacity); _orderBuffer = accelerator.Allocate1D<ulong>(requiredCapacity);
resultBuffer = accelerator.Allocate1D<ulong>(requiredCapacity); _resultBuffer = accelerator.Allocate1D<ulong>(requiredCapacity);
statusBuffer = accelerator.Allocate1D<byte>(requiredCapacity); _statusBuffer = accelerator.Allocate1D<byte>(requiredCapacity);
} }
} }
@@ -496,17 +409,17 @@ public sealed class DivisorCycleCache
var acceleratorIndex = _acceleratorIndex; var acceleratorIndex = _acceleratorIndex;
var stream = AcceleratorStreamPool.Rent(acceleratorIndex); var stream = AcceleratorStreamPool.Rent(acceleratorIndex);
divisorBuffer!.View.CopyFromCPU(stream, divisors); _divisorBuffer!.View.CopyFromCPU(stream, divisors);
powBuffer!.View.CopyFromCPU(stream, powSpan); _powBuffer!.View.CopyFromCPU(stream, powSpan);
orderBuffer!.View.CopyFromCPU(stream, orderSpan); _orderBuffer!.View.CopyFromCPU(stream, orderSpan);
resultBuffer!.View.CopyFromCPU(stream, resultSpan); _resultBuffer!.View.CopyFromCPU(stream, resultSpan);
statusBuffer!.View.CopyFromCPU(stream, statusSpan); _statusBuffer!.View.CopyFromCPU(stream, statusSpan);
do do
{ {
kernel(stream, length, _divisorCyclesBatchSize, divisorBuffer.View, powBuffer.View, orderBuffer.View, resultBuffer.View, statusBuffer.View); kernel(stream, length, _divisorCyclesBatchSize, _divisorBuffer.View, _powBuffer.View, _orderBuffer.View, _resultBuffer.View, _statusBuffer.View);
statusBuffer.View.CopyToCPU(stream, statusSpan); _statusBuffer.View.CopyToCPU(stream, statusSpan);
stream.Synchronize(); stream.Synchronize();
pending = 0; pending = 0;
@@ -520,7 +433,7 @@ public sealed class DivisorCycleCache
} }
while (pending > 0); while (pending > 0);
resultBuffer.View.CopyToCPU(stream, resultSpan); _resultBuffer.View.CopyToCPU(stream, resultSpan);
stream.Synchronize(); stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream); AcceleratorStreamPool.Return(acceleratorIndex, stream);
@@ -7,7 +7,7 @@ namespace PerfectNumbers.Core;
/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path. /// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path.
/// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order. /// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order.
/// </summary> /// </summary>
internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL) public struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL)
{ {
private readonly MontgomeryDivisorData _divisor = divisor; private readonly MontgomeryDivisorData _divisor = divisor;
private readonly ulong _modulus = divisor.Modulus; private readonly ulong _modulus = divisor.Modulus;
@@ -17,10 +17,8 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
private int _gpuBatchSize = GpuConstants.ScanBatchSize; private int _gpuBatchSize = GpuConstants.ScanBatchSize;
// EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards,
// so the synchronization fields from the previous implementation remain commented out here. // so the synchronization fields from the previous implementation remain commented out here.
// private readonly object _sync = new();
private ulong _divisorLimit; private ulong _divisorLimit;
private BigInteger _minK = BigInteger.One; private BigInteger _minK = BigInteger.One;
// private bool _isConfigured;
private readonly ConcurrentFixedCapacityStack<MersenneNumberDivisorByDivisorAccelerator>[] _resourcePool = [.. AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentFixedCapacityStack<MersenneNumberDivisorByDivisorAccelerator>(PerfectNumberConstants.DefaultPoolCapacity))]; private readonly ConcurrentFixedCapacityStack<MersenneNumberDivisorByDivisorAccelerator>[] _resourcePool = [.. AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentFixedCapacityStack<MersenneNumberDivisorByDivisorAccelerator>(PerfectNumberConstants.DefaultPoolCapacity))];
@@ -64,14 +62,8 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
{ {
// EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards,
// so synchronization and runtime configuration guards are unnecessary here. // so synchronization and runtime configuration guards are unnecessary here.
// lock (_sync)
// {
// _divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime);
// _isConfigured = true;
// }
_divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime); _divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime);
// _isConfigured = true;
} }
public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor)
@@ -80,16 +72,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
int batchCapacity; int batchCapacity;
// EvenPerfectBitScanner only calls into this tester after configuring it once, so we can read the cached values without locking. // EvenPerfectBitScanner only calls into this tester after configuring it once, so we can read the cached values without locking.
// lock (_sync)
// {
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
// allowedMax = ComputeAllowedMaxDivisorGpu(prime, _divisorLimit);
// batchCapacity = _gpuBatchSize;
// }
allowedMax = ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); allowedMax = ComputeAllowedMaxDivisorGpu(prime, _divisorLimit);
batchCapacity = _gpuBatchSize; batchCapacity = _gpuBatchSize;
@@ -103,7 +85,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
bool composite; bool composite;
bool coveredRange; bool coveredRange;
ulong processedCount;
ulong lastProcessed; ulong lastProcessed;
// GpuPrimeWorkLimiter(); // GpuPrimeWorkLimiter();
@@ -135,8 +116,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
resources.Counts, resources.Counts,
resources.Cycles, resources.Cycles,
out lastProcessed, out lastProcessed,
out coveredRange, out coveredRange);
out processedCount);
ReturnBatchResources(acceleratorIndex, resources); ReturnBatchResources(acceleratorIndex, resources);
// Monitor.Exit(gpuLease.ExecutionLock); // Monitor.Exit(gpuLease.ExecutionLock);
@@ -213,8 +193,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
int[] counts, int[] counts,
ulong[] cycles, ulong[] cycles,
out ulong lastProcessed, out ulong lastProcessed,
out bool coveredRange, out bool coveredRange)
out ulong processedCount)
{ {
int batchCapacity = (int)divisorDataBuffer.Length; int batchCapacity = (int)divisorDataBuffer.Length;
bool composite = false; bool composite = false;
@@ -244,7 +223,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
if (!hasCandidates) if (!hasCandidates)
{ {
coveredRange = true; coveredRange = true;
processedCount = processedCountLocal;
lastProcessed = lastProcessedLocal; lastProcessed = lastProcessedLocal;
return false; return false;
} }
@@ -272,7 +250,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
ArrayView1D<ulong, Stride1D.Dense> exponentViewDevice = exponentBuffer.View; ArrayView1D<ulong, Stride1D.Dense> exponentViewDevice = exponentBuffer.View;
ArrayView1D<byte, Stride1D.Dense> hitsView = hitsBuffer.View; ArrayView1D<byte, Stride1D.Dense> hitsView = hitsBuffer.View;
ArrayView1D<int, Stride1D.Dense> hitIndexView = hitIndexBuffer.View; ArrayView1D<int, Stride1D.Dense> hitIndexView = hitIndexBuffer.View;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
static bool ProcessRange( static bool ProcessRange(
PrimeOrderCalculatorAccelerator gpu, PrimeOrderCalculatorAccelerator gpu,
ulong prime, ulong prime,
@@ -291,7 +269,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
Span<int> offsetSpan, Span<int> offsetSpan,
Span<int> countSpan, Span<int> countSpan,
Span<ulong> cycleSpan, Span<ulong> cycleSpan,
FixedCapacityStack<MontgomeryDivisorData> divisorPool,
int chunkCountBaseline, int chunkCountBaseline,
int acceleratorIndex, int acceleratorIndex,
ref bool composite, ref bool composite,
@@ -341,8 +318,8 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
if (localStepper.IsAdmissible()) if (localStepper.IsAdmissible())
{ {
ulong divisorValue = (ulong)nextDivisor128; ulong divisorValue = (ulong)nextDivisor128;
MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue); MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
ulong divisorCycle = ResolveDivisorCycle(gpu, divisorValue, prime, in montgomeryData); ulong divisorCycle = ResolveDivisorCycle(gpu, divisorValue, prime, montgomeryData);
if (divisorCycle == prime) if (divisorCycle == prime)
{ {
processedCountLocal += (ulong)(i + 1); processedCountLocal += (ulong)(i + 1);
@@ -351,7 +328,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
return false; return false;
} }
divisorPool.Return(montgomeryData);
divisorSpan[admissibleCount] = divisorValue; divisorSpan[admissibleCount] = divisorValue;
divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue);
offsetSpan[admissibleCount] = admissibleCount; offsetSpan[admissibleCount] = admissibleCount;
@@ -434,7 +410,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
offsetSpan, offsetSpan,
countSpan, countSpan,
cycleSpan, cycleSpan,
divisorPool,
chunkCountBaseline, chunkCountBaseline,
acceleratorIndex, acceleratorIndex,
ref composite, ref composite,
@@ -446,7 +421,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
if (composite) if (composite)
{ {
coveredRange = true; coveredRange = true;
processedCount = processedCountLocal;
lastProcessed = lastProcessedLocal; lastProcessed = lastProcessedLocal;
return true; return true;
} }
@@ -474,7 +448,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
offsetSpan, offsetSpan,
countSpan, countSpan,
cycleSpan, cycleSpan,
divisorPool,
chunkCountBaseline, chunkCountBaseline,
acceleratorIndex, acceleratorIndex,
ref composite, ref composite,
@@ -487,7 +460,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN
processedAll = processedTop && processedBottom; processedAll = processedTop && processedBottom;
coveredRange = composite || processedAll; coveredRange = composite || processedAll;
processedCount = processedCountLocal;
lastProcessed = lastProcessedLocal; lastProcessed = lastProcessedLocal;
return composite; return composite;
} }
@@ -163,7 +163,7 @@ internal static partial class PrimeOrderGpuHeuristics
} }
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
public static GpuPow2ModStatus TryPow2Mod(PrimeOrderCalculatorAccelerator gpu, ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData? divisorData) public static GpuPow2ModStatus TryPow2Mod(PrimeOrderCalculatorAccelerator gpu, ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData divisorData)
{ {
Span<ulong> exponents = stackalloc ulong[1]; Span<ulong> exponents = stackalloc ulong[1];
Span<ulong> remainders = stackalloc ulong[1]; Span<ulong> remainders = stackalloc ulong[1];
@@ -175,7 +175,7 @@ internal static partial class PrimeOrderGpuHeuristics
} }
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
public static GpuPow2ModStatus TryPow2ModBatch(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> exponents, ulong prime, Span<ulong> remainders, in MontgomeryDivisorData? divisorData) public static GpuPow2ModStatus TryPow2ModBatch(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> exponents, ulong prime, Span<ulong> remainders, in MontgomeryDivisorData divisorData)
{ {
// This will never occur in production code // This will never occur in production code
// if (exponents.Length == 0) // if (exponents.Length == 0)
@@ -364,7 +364,7 @@ internal static partial class PrimeOrderGpuHeuristics
return GpuPow2ModStatus.Success; return GpuPow2ModStatus.Success;
} }
private static bool TryComputeOnGpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> exponents, ulong prime, MontgomeryDivisorData? divisorData, Span<ulong> results) private static bool TryComputeOnGpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> exponents, ulong prime, MontgomeryDivisorData divisorData, Span<ulong> results)
{ {
// GpuPrimeWorkLimiter.Acquire(); // GpuPrimeWorkLimiter.Acquire();
int acceleratorIndex = gpu.AcceleratorIndex; int acceleratorIndex = gpu.AcceleratorIndex;
@@ -379,10 +379,11 @@ internal static partial class PrimeOrderGpuHeuristics
var exponentBufferView = gpu.InputView; var exponentBufferView = gpu.InputView;
var remainderBufferView = gpu.OutputUlongView; var remainderBufferView = gpu.OutputUlongView;
if (divisorData is null) // if (divisorData.Equals(MontgomeryDivisorData.Empty) is null)
{ // if (divisorData is null)
divisorData = MontgomeryDivisorData.Empty; // {
} // divisorData = MontgomeryDivisorData.Empty;
// }
AcceleratorStream? stream = AcceleratorStreamPool.Rent(acceleratorIndex); AcceleratorStream? stream = AcceleratorStreamPool.Rent(acceleratorIndex);
exponentBufferView.CopyFromCPU(stream, exponents); exponentBufferView.CopyFromCPU(stream, exponents);
@@ -483,7 +483,7 @@ public sealed class HeuristicCombinedPrimeTester
internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate)
{ {
ulong divisor = candidate.Value; ulong divisor = candidate.Value;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength);
return new(candidate, divisorData, cycleLength, hasCycleHint); return new(candidate, divisorData, cycleLength, hasCycleHint);
} }
+1 -1
View File
@@ -413,7 +413,7 @@ public sealed class HeuristicPrimeTester
internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate)
{ {
ulong divisor = candidate.Value; ulong divisor = candidate.Value;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength);
return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint);
} }
+4 -12
View File
@@ -394,10 +394,8 @@ public class MersenneDivisorCycles
if (divisor <= ulong.MaxValue) if (divisor <= ulong.MaxValue)
{ {
ulong divisor64 = (ulong)divisor; ulong divisor64 = (ulong)divisor;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64);
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor64);
ulong cycle = CalculateCycleLengthCpu(divisor64, divisorData); ulong cycle = CalculateCycleLengthCpu(divisor64, divisorData);
divisorPool.Return(divisorData);
return (UInt128)cycle; return (UInt128)cycle;
} }
@@ -441,10 +439,8 @@ public class MersenneDivisorCycles
if (divisor <= ulong.MaxValue) if (divisor <= ulong.MaxValue)
{ {
ulong divisor64 = (ulong)divisor; ulong divisor64 = (ulong)divisor;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64);
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor64);
ulong cycle = CalculateCycleLengthHybrid(gpu, divisor64, divisorData); ulong cycle = CalculateCycleLengthHybrid(gpu, divisor64, divisorData);
divisorPool.Return(divisorData);
return (UInt128)cycle; return (UInt128)cycle;
} }
@@ -489,10 +485,8 @@ public class MersenneDivisorCycles
if (divisor <= ulong.MaxValue) if (divisor <= ulong.MaxValue)
{ {
ulong divisor64 = (ulong)divisor; ulong divisor64 = (ulong)divisor;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64);
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor64);
ulong cycle = CalculateCycleLengthGpu(gpu, divisor64, divisorData); ulong cycle = CalculateCycleLengthGpu(gpu, divisor64, divisorData);
divisorPool.Return(divisorData);
return (UInt128)cycle; return (UInt128)cycle;
} }
@@ -1168,7 +1162,6 @@ public class MersenneDivisorCycles
// (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength)); // (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength));
int localCycleIndex = 0; int localCycleIndex = 0;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
for (ulong divisor = localStart; divisor <= localEnd; divisor++) for (ulong divisor = localStart; divisor <= localEnd; divisor++)
{ {
@@ -1188,9 +1181,8 @@ public class MersenneDivisorCycles
// TODO: Migrate this generation path to the shared // TODO: Migrate this generation path to the shared
// unrolled-hex calculator so the CPU generator matches // unrolled-hex calculator so the CPU generator matches
// the GPU benchmark leader for large divisors. // the GPU benchmark leader for large divisors.
MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
ulong cycle = CalculateCycleLengthGpu(gpu, divisor, divisorData); ulong cycle = CalculateCycleLengthGpu(gpu, divisor, divisorData);
divisorPool.Return(divisorData);
localCycles[localCycleIndex++] = (divisor, cycle); localCycles[localCycleIndex++] = (divisor, cycle);
} }
+35 -13
View File
@@ -1,24 +1,46 @@
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core; namespace PerfectNumbers.Core;
public sealed class MontgomeryDivisorData public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared)
{ {
public static readonly MontgomeryDivisorData Empty = new(); public static readonly MontgomeryDivisorData Empty = new();
public ulong Modulus;
public ulong NPrime;
public ulong MontgomeryOne;
public ulong MontgomeryTwo;
public ulong MontgomeryTwoSquared;
public MontgomeryDivisorData() [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static ulong ComputeMontgomeryResidue(UInt128 value, ulong modulus) => (ulong)(value % modulus);
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{ {
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
} }
public MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared) [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static MontgomeryDivisorData FromModulus(ulong modulus)
{ {
Modulus = modulus; ulong nPrime = ComputeMontgomeryNPrime(modulus);
NPrime = nPrime; ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
MontgomeryOne = montgomeryOne; ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus);
MontgomeryTwo = montgomeryTwo; ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
MontgomeryTwoSquared = montgomeryTwoSquared;
return new(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
} }
public readonly ulong Modulus = modulus;
public readonly ulong NPrime = nPrime;
public readonly ulong MontgomeryOne = montgomeryOne;
public readonly ulong MontgomeryTwo = montgomeryTwo;
public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared;
} }
@@ -1,86 +0,0 @@
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
public static class MontgomeryDivisorDataPool
{
[ThreadStatic]
private static FixedCapacityStack<MontgomeryDivisorData>? _pool;
public static FixedCapacityStack<MontgomeryDivisorData> Shared
{
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
get => _pool ??= new(PerfectNumberConstants.DefaultPoolCapacity);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static ulong ComputeMontgomeryResidue(UInt128 value, ulong modulus) => (ulong)(value % modulus);
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static MontgomeryDivisorData Rent(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared)
{
var pool = _pool ??= new(PerfectNumberConstants.DefaultPoolCapacity);
if (pool.Pop() is { } pooled)
{
pooled.Modulus = modulus;
pooled.NPrime = nPrime;
pooled.MontgomeryOne = montgomeryOne;
pooled.MontgomeryTwo = montgomeryTwo;
pooled.MontgomeryTwoSquared = montgomeryTwoSquared;
return pooled;
}
return new(modulus, nPrime, montgomeryOne, montgomeryTwo, montgomeryTwoSquared);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static void Return(MontgomeryDivisorData divisorData) => _pool!.Push(divisorData);
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static MontgomeryDivisorData FromModulus(this FixedCapacityStack<MontgomeryDivisorData> queue, ulong modulus)
{
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return queue.Rent(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static MontgomeryDivisorData Rent(this FixedCapacityStack<MontgomeryDivisorData> queue, ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared)
{
if (queue.Pop() is { } pooled)
{
pooled.Modulus = modulus;
pooled.NPrime = nPrime;
pooled.MontgomeryOne = montgomeryOne;
pooled.MontgomeryTwo = montgomeryTwo;
pooled.MontgomeryTwoSquared = montgomeryTwoSquared;
return pooled;
}
return new(modulus, nPrime, montgomeryOne, montgomeryTwo, montgomeryTwoSquared);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static void Return(this FixedCapacityStack<MontgomeryDivisorData> queue, MontgomeryDivisorData divisorData) => queue.Push(divisorData);
}
+9 -8
View File
@@ -2,7 +2,7 @@ using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core namespace PerfectNumbers.Core
{ {
public sealed class PartialFactorResult public sealed class PartialFactorResult(in MontgomeryDivisorData divisorData)
{ {
[ThreadStatic] [ThreadStatic]
private static PartialFactorResult? s_poolHead; private static PartialFactorResult? s_poolHead;
@@ -13,26 +13,27 @@ namespace PerfectNumbers.Core
public readonly List<PartialFactorPendingEntry> PendingFactors = new (PrimeOrderConstants.GpuSmallPrimeFactorSlots); public readonly List<PartialFactorPendingEntry> PendingFactors = new (PrimeOrderConstants.GpuSmallPrimeFactorSlots);
public readonly ulong[] StackCandidates = new ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots << 3]; public readonly ulong[] StackCandidates = new ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots << 3];
public readonly bool[] StackEvaluations = new bool[PrimeOrderConstants.GpuSmallPrimeFactorSlots << 3]; public readonly bool[] StackEvaluations = new bool[PrimeOrderConstants.GpuSmallPrimeFactorSlots << 3];
public ExponentRemainderStepperCpu ExponentRemainderStepper = new(divisorData);
private PartialFactorResult? _next; private PartialFactorResult? _next;
public ulong Cofactor; public ulong Cofactor;
public bool FullyFactored; public bool FullyFactored;
public int Count; public int Count;
public bool CofactorIsPrime; public bool CofactorIsPrime;
private PartialFactorResult() public static PartialFactorResult Rent(in MontgomeryDivisorData divisorData)
{
}
public static PartialFactorResult Rent()
{ {
PartialFactorResult? instance = s_poolHead; PartialFactorResult? instance = s_poolHead;
if (instance is null) if (instance is null)
{ {
instance = new PartialFactorResult(); instance = new PartialFactorResult(divisorData);
} }
else else
{ {
s_poolHead = instance._next; s_poolHead = instance._next;
if (!instance.ExponentRemainderStepper.MatchesDivisor(divisorData, 0UL))
{
instance.ExponentRemainderStepper = new(divisorData, 0UL);
}
} }
instance._next = null; instance._next = null;
@@ -43,7 +44,7 @@ namespace PerfectNumbers.Core
return instance; return instance;
} }
public static readonly PartialFactorResult Empty = new() public static readonly PartialFactorResult Empty = new(MontgomeryDivisorData.Empty)
{ {
Cofactor = 1UL, Cofactor = 1UL,
FullyFactored = true FullyFactored = true
@@ -25,7 +25,7 @@ internal static partial class PrimeOrderCalculator
ulong phi = prime - 1UL; ulong phi = prime - 1UL;
PartialFactorResult phiFactors = PartialFactorCpu(phi, config); PartialFactorResult phiFactors = PartialFactorCpu(phi, divisorData, config);
ulong result; ulong result;
if (phiFactors.Factors is null) if (phiFactors.Factors is null)
@@ -46,7 +46,6 @@ internal static partial class PrimeOrderCalculator
in UInt128? previousOrder, in UInt128? previousOrder,
in PrimeOrderCalculatorConfig config) in PrimeOrderCalculatorConfig config)
{ {
MontgomeryDivisorData divisorData;
UInt128 result; UInt128 result;
if (prime <= ulong.MaxValue) if (prime <= ulong.MaxValue)
{ {
@@ -65,16 +64,13 @@ internal static partial class PrimeOrderCalculator
} }
ulong prime64 = (ulong)prime; ulong prime64 = (ulong)prime;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime64);
divisorData = divisorPool.FromModulus(prime64);
ulong order64 = CalculateCpu(prime64, previous, divisorData, config); ulong order64 = CalculateCpu(prime64, previous, divisorData, config);
divisorPool.Return(divisorData);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
} }
else else
{ {
divisorData = MontgomeryDivisorData.Empty; result = CalculateWideInternalCpu(prime, previousOrder, MontgomeryDivisorData.Empty, config);
result = CalculateWideInternalCpu(prime, previousOrder, divisorData, config);
} }
return result; return result;
@@ -561,7 +557,7 @@ internal static partial class PrimeOrderCalculator
SortFactorExponentSpans(factorsSpan, exponentsSpan, factorCount); SortFactorExponentSpans(factorsSpan, exponentsSpan, factorCount);
ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); ExponentRemainderStepperCpu stepper = factors.ExponentRemainderStepper;
Span<ulong> stackCandidates = new(factors.StackCandidates); Span<ulong> stackCandidates = new(factors.StackCandidates);
Span<bool> stackEvaluations = new(factors.StackEvaluations); Span<bool> stackEvaluations = new(factors.StackEvaluations);
@@ -579,7 +575,6 @@ internal static partial class PrimeOrderCalculator
ProcessExponentLoweringPrime(stackCandidates[..exponent], stackEvaluations[..exponent], ref order, primeFactor, exponent, ref stepper); ProcessExponentLoweringPrime(stackCandidates[..exponent], stackEvaluations[..exponent], ref order, primeFactor, exponent, ref stepper);
} }
ThreadStaticPools.ReturnExponentStepperCpu(stepper);
return order; return order;
} }
@@ -673,7 +668,7 @@ internal static partial class PrimeOrderCalculator
return false; return false;
} }
PartialFactorResult factorization = PartialFactorCpu(order, config); PartialFactorResult factorization = PartialFactorCpu(order, divisorData, config);
if (factorization.Factors is null) if (factorization.Factors is null)
{ {
factorization.Dispose(); factorization.Dispose();
@@ -810,7 +805,7 @@ internal static partial class PrimeOrderCalculator
} }
// Reuse the partial factorization from TryConfirmOrderCpu when available. // Reuse the partial factorization from TryConfirmOrderCpu when available.
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorCpu(order, config); PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorCpu(order, divisorData, config);
try try
{ {
if (orderFactors.Factors is null) if (orderFactors.Factors is null)
@@ -1106,7 +1101,7 @@ internal static partial class PrimeOrderCalculator
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget)
{ {
PartialFactorResult factorization = PartialFactorCpu(candidate, config); PartialFactorResult factorization = PartialFactorCpu(candidate, divisorData, config);
if (factorization.Factors is null) if (factorization.Factors is null)
{ {
@@ -1270,7 +1265,7 @@ internal static partial class PrimeOrderCalculator
=> exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; => exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL;
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static PartialFactorResult PartialFactorCpu(ulong value, in PrimeOrderCalculatorConfig config) private static PartialFactorResult PartialFactorCpu(ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config)
{ {
if (value <= 1UL) if (value <= 1UL)
{ {
@@ -1279,7 +1274,7 @@ internal static partial class PrimeOrderCalculator
// stackalloc is faster than pooling // stackalloc is faster than pooling
var result = PartialFactorResult.Rent(); var result = PartialFactorResult.Rent(divisorData);
ulong[] factors = result.Factors; ulong[] factors = result.Factors;
Span<ulong> primeSlots = new(factors); Span<ulong> primeSlots = new(factors);
int[] exponents = result.Exponents; int[] exponents = result.Exponents;
+15 -17
View File
@@ -76,7 +76,7 @@ internal static partial class PrimeOrderCalculator
return gpuOrder; return gpuOrder;
} }
PartialFactorResult phiFactors = PartialFactorGpu(gpu, phi, config); PartialFactorResult phiFactors = PartialFactorGpu(gpu, phi, divisorData, config);
ulong result; ulong result;
if (phiFactors.Factors is null) if (phiFactors.Factors is null)
@@ -117,10 +117,8 @@ internal static partial class PrimeOrderCalculator
} }
ulong prime64 = (ulong)prime; ulong prime64 = (ulong)prime;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; divisorData = MontgomeryDivisorData.FromModulus(prime64);
divisorData = divisorPool.FromModulus(prime64);
ulong order64 = CalculateGpu(gpu, prime64, previous, divisorData, config); ulong order64 = CalculateGpu(gpu, prime64, previous, divisorData, config);
divisorPool.Return(divisorData);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
} }
else else
@@ -133,7 +131,7 @@ internal static partial class PrimeOrderCalculator
} }
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static UInt128 CalculateWideInternalGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) private static UInt128 CalculateWideInternalGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config)
{ {
if (prime <= UInt128.One) if (prime <= UInt128.One)
{ {
@@ -233,7 +231,7 @@ internal static partial class PrimeOrderCalculator
return result != 0; return result != 0;
} }
private static UInt128 ExponentLoweringWideGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) private static UInt128 ExponentLoweringWideGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
{ {
ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool; ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool;
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors; ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
@@ -281,7 +279,7 @@ internal static partial class PrimeOrderCalculator
return order; return order;
} }
private static UInt128 FinishStrictlyWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData? divisorData) private static UInt128 FinishStrictlyWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData divisorData)
{ {
UInt128 phi = prime - UInt128.One; UInt128 phi = prime - UInt128.One;
Dictionary<UInt128, int> counts = new(capacity: 8); Dictionary<UInt128, int> counts = new(capacity: 8);
@@ -320,7 +318,7 @@ internal static partial class PrimeOrderCalculator
return order; return order;
} }
private static UInt128 InitializeStartingOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData? divisorData) private static UInt128 InitializeStartingOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData divisorData)
{ {
UInt128 order = phi; UInt128 order = phi;
UInt128 mod8 = prime & (UInt128)7UL; UInt128 mod8 = prime & (UInt128)7UL;
@@ -337,7 +335,7 @@ internal static partial class PrimeOrderCalculator
} }
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config)
{ {
if (value <= 1UL) if (value <= 1UL)
{ {
@@ -345,7 +343,7 @@ internal static partial class PrimeOrderCalculator
} }
// stackalloc is faster than pooling // stackalloc is faster than pooling
var result = PartialFactorResult.Rent(); var result = PartialFactorResult.Rent(divisorData);
ulong[] factors = result.Factors; ulong[] factors = result.Factors;
Span<ulong> primeSlots = new (factors); Span<ulong> primeSlots = new (factors);
int[] exponents = result.Exponents; int[] exponents = result.Exponents;
@@ -478,7 +476,7 @@ internal static partial class PrimeOrderCalculator
} }
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static UInt128 Pow2ModWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData? divisorData) private static UInt128 Pow2ModWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData)
{ {
if (modulus == UInt128.One) if (modulus == UInt128.One)
{ {
@@ -509,7 +507,7 @@ internal static partial class PrimeOrderCalculator
PrimeOrderCalculatorAccelerator gpu, PrimeOrderCalculatorAccelerator gpu,
in UInt128 prime, in UInt128 prime,
in UInt128? previousOrder, in UInt128? previousOrder,
in MontgomeryDivisorData? divisorData, in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config, in PrimeOrderCalculatorConfig config,
in UInt128 phi, in UInt128 phi,
in PartialFactorResult128 phiFactors) in PartialFactorResult128 phiFactors)
@@ -540,7 +538,7 @@ internal static partial class PrimeOrderCalculator
return candidateOrder; return candidateOrder;
} }
private static bool TryConfirmCandidateWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData? divisorData) private static bool TryConfirmCandidateWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData)
{ {
PartialFactorResult128 factorization = PartialFactorWide(candidate, config); PartialFactorResult128 factorization = PartialFactorWide(candidate, config);
if (factorization.Factors is null) if (factorization.Factors is null)
@@ -593,7 +591,7 @@ internal static partial class PrimeOrderCalculator
return true; return true;
} }
private static bool TryConfirmOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) private static bool TryConfirmOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config)
{ {
if (order == UInt128.Zero) if (order == UInt128.Zero)
{ {
@@ -670,7 +668,7 @@ internal static partial class PrimeOrderCalculator
} }
// Reuse the partial factorization from TryConfirmOrderCpu when available. // Reuse the partial factorization from TryConfirmOrderCpu when available.
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, config); PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, divisorData, config);
try try
{ {
if (orderFactors.Factors is null) if (orderFactors.Factors is null)
@@ -854,7 +852,7 @@ internal static partial class PrimeOrderCalculator
in UInt128 prime, in UInt128 prime,
in UInt128 order, in UInt128 order,
in UInt128? previousOrder, in UInt128? previousOrder,
in MontgomeryDivisorData? divisorData, in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config, in PrimeOrderCalculatorConfig config,
out UInt128 result) out UInt128 result)
{ {
@@ -1023,7 +1021,7 @@ internal static partial class PrimeOrderCalculator
return true; return true;
} }
private static bool TrySpecialMaxWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) private static bool TrySpecialMaxWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
{ {
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors; ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
int length = factors.Count; int length = factors.Count;
@@ -25,7 +25,7 @@ internal static partial class PrimeOrderCalculator
ulong phi = prime - 1UL; ulong phi = prime - 1UL;
PartialFactorResult phiFactors = PartialFactorHybrid(gpu, phi, config); PartialFactorResult phiFactors = PartialFactorHybrid(gpu, phi, divisorData, config);
ulong result; ulong result;
if (phiFactors.Factors is null) if (phiFactors.Factors is null)
@@ -67,10 +67,8 @@ internal static partial class PrimeOrderCalculator
} }
ulong prime64 = (ulong)prime; ulong prime64 = (ulong)prime;
FixedCapacityStack<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared; divisorData = MontgomeryDivisorData.FromModulus(prime64);
divisorData = divisorPool.FromModulus(prime64);
ulong order64 = CalculateHybrid(gpu, prime64, previous, divisorData, config); ulong order64 = CalculateHybrid(gpu, prime64, previous, divisorData, config);
divisorPool.Return(divisorData);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
} }
else else
@@ -134,7 +132,7 @@ internal static partial class PrimeOrderCalculator
} }
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config)
{ {
if (value <= 1UL) if (value <= 1UL)
{ {
@@ -143,7 +141,7 @@ internal static partial class PrimeOrderCalculator
// stackalloc is faster than pooling // stackalloc is faster than pooling
var result = PartialFactorResult.Rent(); var result = PartialFactorResult.Rent(divisorData);
ulong[] factors = result.Factors; ulong[] factors = result.Factors;
Span<ulong> primeSlots = new(factors); Span<ulong> primeSlots = new(factors);
int[] exponents = result.Exponents; int[] exponents = result.Exponents;
@@ -369,7 +367,7 @@ internal static partial class PrimeOrderCalculator
[MethodImpl(MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget)
{ {
PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, config); PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, divisorData, config);
if (factorization.Factors is null) if (factorization.Factors is null)
{ {
@@ -466,7 +464,7 @@ internal static partial class PrimeOrderCalculator
return false; return false;
} }
PartialFactorResult factorization = PartialFactorHybrid(gpu, order, config); PartialFactorResult factorization = PartialFactorHybrid(gpu, order, divisorData, config);
if (factorization.Factors is null) if (factorization.Factors is null)
{ {
factorization.Dispose(); factorization.Dispose();
@@ -514,7 +512,7 @@ internal static partial class PrimeOrderCalculator
} }
// Reuse the partial factorization from TryConfirmOrderCpu when available. // Reuse the partial factorization from TryConfirmOrderCpu when available.
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, config); PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, divisorData, config);
try try
{ {
if (orderFactors.Factors is null) if (orderFactors.Factors is null)
-9
View File
@@ -182,15 +182,6 @@ namespace PerfectNumbers.Core
get => _bytePool ??= ArrayPool<byte>.Create(); get => _bytePool ??= ArrayPool<byte>.Create();
} }
[ThreadStatic]
private static ArrayPool<MontgomeryDivisorData>? _montgomeryDivisorDataPool;
public static ArrayPool<MontgomeryDivisorData> MontgomeryDivisorDataPool
{
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
get => _montgomeryDivisorDataPool ??= ArrayPool<MontgomeryDivisorData>.Create();
}
[ThreadStatic] [ThreadStatic]
private static ArrayPool<GpuDivisorPartialData>? _gpuDivisorPartialDataPool; private static ArrayPool<GpuDivisorPartialData>? _gpuDivisorPartialDataPool;
+20 -20
View File
@@ -510,38 +510,38 @@ public static class UInt128Extensions
ulong rightLow = (ulong)right; ulong rightLow = (ulong)right;
ulong rightHigh = (ulong)(right >> 64); ulong rightHigh = (ulong)(right >> 64);
var partial0 = Multiply64(leftLow, rightLow); (ulong carry, p0) = Multiply64(leftLow, rightLow);
var partial1 = Multiply64(leftLow, rightHigh); (ulong partial1High, ulong partial1Low) = Multiply64(leftLow, rightHigh);
var partial2 = Multiply64(leftHigh, rightLow);
var partial3 = Multiply64(leftHigh, rightHigh);
p0 = partial0.Low; ulong sum = carry + partial1Low;
ulong carry = partial0.High; carry = sum < partial1Low ? 1UL : 0UL;
ulong sum = carry + partial1.Low;
ulong carryMid = sum < partial1.Low ? 1UL : 0UL; (ulong partial2High, ulong partial2Low) = Multiply64(leftHigh, rightLow);
sum += partial2.Low; sum += partial2Low;
if (sum < partial2.Low) if (sum < partial2Low)
{ {
carryMid++; carry++;
} }
p1 = sum; p1 = sum;
sum = partial1.High + partial2.High; sum = partial1High + partial2High;
ulong carryHigh = sum < partial2.High ? 1UL : 0UL; partial2High = sum < partial2High ? 1UL : 0UL;
sum += partial3.Low;
if (sum < partial3.Low) (partial1High, partial1Low) = Multiply64(leftHigh, rightHigh);
sum += partial1Low;
if (sum < partial1Low)
{ {
carryHigh++; partial2High++;
} }
sum += carryMid; sum += carry;
if (sum < carryMid) if (sum < carry)
{ {
carryHigh++; partial2High++;
} }
p2 = sum; p2 = sum;
p3 = partial3.High + carryHigh; p3 = partial1High + partial2High;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]