From 070e2a14fdefe7a05992e40dd36da6f15bdb062f Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Mon, 24 Nov 2025 15:14:31 +0100 Subject: [PATCH] MontgomeryDivisorData converted to a pooled class to limit GC pressure. MontgomeryDivisorDataGpu added for GPU operations. PrimeOrderCalculator kernels precompiled during warmup. --- .../Pow2ModBenchmarks.cs | 6 +- .../Pow2MontgomeryModBenchmarks.cs | 2 +- .../Pow2MontgomeryModWindowedBenchmarkData.cs | 2 +- .../Gpu/PrimeOrderGpuHeuristicsTests.cs | 18 ++--- .../Gpu/ULongExtensionsGpuTests.cs | 4 +- .../MersenneDivisorCyclesGenerateGpuTests.cs | 4 +- .../MersenneDivisorCyclesTests.cs | 2 +- .../MersenneNumberDivisorCpuTesterTests.cs | 12 +-- .../MersenneNumberDivisorGpuTesterTests.cs | 24 +++--- .../PrimeOrderCalculatorTests.cs | 12 +-- .../Cpu/MersenneCpuDivisorScanSession.cs | 8 +- ...MersenneNumberDivisorByDivisorCpuTester.cs | 58 ++++++-------- PerfectNumbers.Core/DivisorCycleCache.cs | 18 +++-- .../PrimeOrderCalculatorAccelerator.cs | 63 +++++++++------ .../Gpu/Kernels/Pow2MontgomeryKernels.cs | 14 +++- .../Gpu/Kernels/PrimeOrderKernels.cs | 32 +++++--- .../Kernels/PrimeOrderSpecialMaxKernels.cs | 4 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 5 +- .../Gpu/MontgomeryDivisorDataGpu.cs | 10 +++ .../Gpu/Pow2MontgomeryGpuCalculator.cs | 6 +- .../Gpu/PrimeOrderGpuHeuristics.cs | 23 +++--- .../HeuristicCombinedPrimeTester.cs | 2 +- PerfectNumbers.Core/HeuristicPrimeTester.cs | 2 +- PerfectNumbers.Core/MersenneDivisorCycles.cs | 8 +- PerfectNumbers.Core/MontgomeryDivisorData.cs | 59 +++++--------- .../MontgomeryDivisorDataPool.cs | 80 +++++++++++++++++++ PerfectNumbers.Core/PerfectNumberConstants.cs | 2 +- .../PrimeOrderCalculator.Cpu.cs | 2 +- .../PrimeOrderCalculator.Gpu.cs | 4 +- .../PrimeOrderCalculator.Shared.cs | 2 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 6 +- PerfectNumbers.Core/ULongExtensions.Gpu.cs | 8 +- 32 files changed, 310 insertions(+), 192 deletions(-) create mode 100644 PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs create mode 100644 PerfectNumbers.Core/MontgomeryDivisorDataPool.cs diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs index 6f51747..ea0a4c7 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs @@ -470,7 +470,7 @@ public class Pow2ModBenchmarks Console.WriteLine($"Calculating small cycle {smallModulus}"); smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); smallExponents[i] = NextSmallExponent(random); - smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus, MontgomeryDivisorData.FromModulus(smallModulus)); + smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus, MontgomeryDivisorDataPool.Shared.FromModulus(smallModulus)); (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(random); Console.WriteLine($"Calculating large cycle {largeModulus}"); @@ -582,7 +582,7 @@ public class Pow2ModBenchmarks private static ulong CalculateCycleLengthWithHeuristics(ulong modulus, ref ulong? previousPrimeOrder) { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); if (modulus <= 1UL || (modulus & 1UL) == 0UL) { return MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData); @@ -651,7 +651,7 @@ public class Pow2ModBenchmarks private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL) { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); if (modulus <= 1UL || (modulus & 1UL) == 0UL) { return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData); diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index b9323ef..5d11eda 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -37,7 +37,7 @@ public class Pow2MontgomeryModBenchmarks ulong smallModulus = NextSmallOddModulus(); _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); _smallExponents[i] = NextSmallExponent(); - _smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus, MontgomeryDivisorData.FromModulus(smallModulus)); + _smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus, MontgomeryDivisorDataPool.Shared.FromModulus(smallModulus)); Console.WriteLine($"Calculating large cycle {i + 1} with heuristics"); (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(); diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs index 1f01af9..37a1c18 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs @@ -64,7 +64,7 @@ internal static class Pow2MontgomeryModWindowedBenchmarkInputsProvider { Pow2MontgomeryModWindowedBenchmarkSeed seed = seeds[i]; ulong modulus = 2UL * seed.Multiplier * seed.Exponent + 1UL; - MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + MontgomeryDivisorData divisor = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); ulong cycleLength = MersenneDivisorCycles.CalculateCycleLength(modulus, divisor); ulong reducedExponent = seed.Exponent % cycleLength; cases[i] = new Pow2MontgomeryModWindowedBenchmarkCase( diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs index f9da5bb..5ac6ec4 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs @@ -23,7 +23,7 @@ public class PrimeOrderGpuHeuristicsTests try { PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); @@ -40,7 +40,7 @@ public class PrimeOrderGpuHeuristicsTests { const ulong exponent = 138_000_001UL; const ulong prime = 138_000_001UL; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Success); @@ -55,7 +55,7 @@ public class PrimeOrderGpuHeuristicsTests PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); try { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); @@ -76,7 +76,7 @@ public class PrimeOrderGpuHeuristicsTests try { const ulong prime = 193UL; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(16UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); @@ -102,7 +102,7 @@ public class PrimeOrderGpuHeuristicsTests ulong[] remainders = new ulong[2]; ulong[] exponents = { 1UL, 2UL }; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); @@ -121,7 +121,7 @@ public class PrimeOrderGpuHeuristicsTests ulong[] remainders = { 123UL, 456UL }; ulong[] exponents = { 138_000_001UL, 138_000_011UL }; const ulong prime = 138_000_001UL; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); ulong expectedFirst = (ulong)BigInteger.ModPow(2, exponents[0], prime); @@ -145,7 +145,7 @@ public class PrimeOrderGpuHeuristicsTests ulong[] exponents = { 1UL, 2UL }; ulong[] remainders = { 5UL, 6UL }; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); @@ -170,7 +170,7 @@ public class PrimeOrderGpuHeuristicsTests ulong[] remainders = { 9UL, 11UL }; const ulong prime = 193UL; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); @@ -362,7 +362,7 @@ public class PrimeOrderGpuHeuristicsTests { PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(prime); ulong result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); result.Should().Be(3UL); diff --git a/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs index 789ab31..89c6e38 100644 --- a/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs @@ -17,7 +17,7 @@ public class ULongExtensionsGpuTests [InlineData(65535UL, 2147483647UL)] public void Pow2ModWindowedGpu_matches_montgomery_cpu_for_known_values(ulong exponent, ulong modulus) { - MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + MontgomeryDivisorData divisor = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); ulong expected = exponent.Pow2MontgomeryModWindowedCpu(divisor, keepMontgomery: false); exponent.Pow2ModWindowedGpu(modulus).Should().Be(expected); @@ -51,7 +51,7 @@ public class ULongExtensionsGpuTests foreach (ulong modulus in moduli) { - MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + MontgomeryDivisorData divisor = MontgomeryDivisorDataPool.Shared.FromModulus(modulus); foreach (ulong exponent in exponents) { ulong expected = exponent.Pow2MontgomeryModWindowedCpu(divisor, keepMontgomery: false); diff --git a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs index af63152..58a2967 100644 --- a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs @@ -38,7 +38,7 @@ public class MersenneDivisorCyclesGenerateGpuTests foreach (var pair in pairs) { expected.Should().Contain(pair.divisor); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(pair.divisor); pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData)); } } @@ -87,7 +87,7 @@ public class MersenneDivisorCyclesGenerateGpuTests foreach (var pair in pairs) { uniqueDivisors.Add(pair.divisor).Should().BeTrue(); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(pair.divisor); pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData)); } diff --git a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs index 80c844a..a332fdc 100644 --- a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs @@ -18,7 +18,7 @@ public class MersenneDivisorCyclesTests [MemberData(nameof(Divisors))] public void CalculateCycleLength_returns_expected_value(ulong divisor, ulong expected) { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(divisor); MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData).Should().Be(expected); MersenneDivisorCycles.CalculateCycleLengthGpu(divisor).Should().Be(expected); } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index 87ff295..5581905 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -39,13 +39,13 @@ public class MersenneNumberDivisorCpuTesterTests ulong[] primes = { 5UL, 7UL, 11UL, 13UL }; byte[] hits = new byte[primes.Length]; - ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL)); - session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits); + ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorDataPool.Shared.FromModulus(23UL)); + session.CheckDivisor(23UL, MontgomeryDivisorDataPool.Shared.FromModulus(23UL), cycle23, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); Array.Fill(hits, (byte)0); - ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL)); - session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits); + ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorDataPool.Shared.FromModulus(31UL)); + session.CheckDivisor(31UL, MontgomeryDivisorDataPool.Shared.FromModulus(31UL), cycle31, primes, hits); hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 }); } finally @@ -67,7 +67,7 @@ public class MersenneNumberDivisorCpuTesterTests ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(7UL); ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData); session.CheckDivisor(7UL, divisorData, cycle, exponents, hits); @@ -92,7 +92,7 @@ public class MersenneNumberDivisorCpuTesterTests ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(11UL); ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData); session.CheckDivisor(11UL, divisorData, cycle, exponents, hits); diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 707ec22..8cb8323 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -170,7 +170,7 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL }; byte[] hits = new byte[primes.Length]; - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes, hits); + session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), 0UL, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0 }); } @@ -193,13 +193,13 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL, 43UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 }); Array.Fill(hits, (byte)0); - ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL)); - session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits); + ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorDataPool.Shared.FromModulus(13367UL)); + session.CheckDivisor(13367UL, MontgomeryDivisorDataPool.Shared.FromModulus(13367UL), cycle13367, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); } finally @@ -225,8 +225,8 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 }); } @@ -252,7 +252,7 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(7UL); ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData); session.CheckDivisor(7UL, divisorData, cycle, exponents, hits); @@ -281,7 +281,7 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL }; byte[] hits = new byte[exponents.Length]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(11UL); ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData); session.CheckDivisor(11UL, divisorData, cycle, exponents, hits); @@ -310,7 +310,7 @@ public class MersenneNumberDivisorGpuTesterTests var patchedTable = new List<(ulong Divisor, ulong Cycle)>(); ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(191UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(191UL); patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLength(191UL, divisorData); tableField.SetValue(cycles, patchedTable); @@ -353,8 +353,8 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorDataPool.Shared.FromModulus(223UL), cycle223, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 }); } diff --git a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs index f3a27d3..19eaaa8 100644 --- a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs @@ -46,7 +46,7 @@ public class PrimeOrderCalculatorTests ulong result = PrimeOrderCalculator.Calculate( prime, previousOrder: null, - MontgomeryDivisorData.FromModulus(prime), + MontgomeryDivisorDataPool.Shared.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); @@ -63,14 +63,14 @@ public class PrimeOrderCalculatorTests ulong heuristic = PrimeOrderCalculator.Calculate( prime, previousOrder, - MontgomeryDivisorData.FromModulus(prime), + MontgomeryDivisorDataPool.Shared.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); ulong strict = PrimeOrderCalculator.Calculate( prime, previousOrder: null, - MontgomeryDivisorData.FromModulus(prime), + MontgomeryDivisorDataPool.Shared.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); @@ -86,7 +86,7 @@ public class PrimeOrderCalculatorTests ulong result = PrimeOrderCalculator.Calculate( prime, previousOrder: null, - MontgomeryDivisorData.FromModulus(prime), + MontgomeryDivisorDataPool.Shared.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); @@ -106,7 +106,7 @@ public class PrimeOrderCalculatorTests ulong result = PrimeOrderCalculator.Calculate( prime: 13UL, previousOrder: null, - MontgomeryDivisorData.FromModulus(13UL), + MontgomeryDivisorDataPool.Shared.FromModulus(13UL), config, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu); @@ -123,7 +123,7 @@ public class PrimeOrderCalculatorTests maxPowChecks: 1, mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(239UL); ulong heuristic = PrimeOrderCalculator.Calculate( prime: 239UL, previousOrder: null, diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs index 44984ea..2d36028 100644 --- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs +++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs @@ -21,10 +21,10 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi // return; // } - MontgomeryDivisorData cachedData = divisorData; - if (cachedData.Modulus != divisor) + if (divisorData.Modulus != divisor) { - cachedData = MontgomeryDivisorData.FromModulus(divisor); + throw new InvalidOperationException("Divisor data is for a different modulus"); + // cachedData = MontgomeryDivisorData.FromModulus(divisor); } if (divisorCycle == 0UL) @@ -39,7 +39,7 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi // Keep these remainder steppers in place so future updates continue reusing the previously computed residues. // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. - var exponentStepper = new ExponentRemainderStepperCpu(cachedData); + var exponentStepper = new ExponentRemainderStepperCpu(divisorData); var cycleStepper = new CycleRemainderStepper(divisorCycle); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index e497ee1..5b3ba63 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -204,43 +204,37 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - + var divisorPool = MontgomeryDivisorDataPool.Shared; while (divisor.CompareTo(limit) <= 0) { bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0; if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) { ulong candidate = divisor.Low; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + MontgomeryDivisorData divisorData = divisorPool.FromModulus(candidate); ulong divisorCycle; // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. - // if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) - // { - // divisorCycle = cycleCache.GetCycleLength(candidate); - // } - // else + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + gpu, + candidate, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) { - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( - gpu, - candidate, - prime, - divisorData, - out ulong computedCycle, - out bool primeOrderFailed) || computedCycle == 0UL) - { - // Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles - // for the exponents scanned here, so skip the unused cache fallback and compute directly. - divisorCycle = MersenneDivisorCycles.CalculateCycleLength( - candidate, - divisorData, - skipPrimeOrderHeuristic: primeOrderFailed); - } - else - { - divisorCycle = computedCycle; - } + // Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles + // for the exponents scanned here, so skip the unused cache fallback and compute directly. + divisorCycle = MersenneDivisorCycles.CalculateCycleLength( + candidate, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; } + divisorPool.Return(divisorData); if (divisorCycle == prime) { // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides @@ -290,18 +284,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv out bool processedAll, out ulong foundDivisor) { - processedAll = false; - foundDivisor = 0UL; - bool canAdvance = step <= limit; var gpu = PrimeOrderCalculatorAccelerator.Rent(1); + Queue divisorPool = MontgomeryDivisorDataPool.Shared; if (!canAdvance) { if (divisor <= limit) { if (remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); ulong divisorCycle; if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( gpu, @@ -321,6 +313,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv divisorCycle = computedCycle; } + divisorPool.Return(divisorData); + if (divisorCycle == prime) { foundDivisor = divisor; @@ -351,13 +345,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } ulong remainingIterations = ((limit - divisor) / step) + 1UL; - while (true) { bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0; if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); ulong divisorCycle; if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( gpu, @@ -377,6 +370,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv divisorCycle = computedCycle; } + divisorPool.Return(divisorData); if (divisorCycle == prime) { foundDivisor = divisor; diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 39c0d39..3f8f571 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -76,19 +76,23 @@ public sealed class DivisorCycleCache // } ulong[] snapshot = _snapshot; + ulong cycleLength; if (divisor < (ulong)snapshot.Length) { - ulong cached = snapshot[divisor]; - if (cached == 0UL) + cycleLength = snapshot[divisor]; + if (cycleLength == 0UL) { throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } - return cached; + return cycleLength; } - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: true); + Queue divisorPool = MontgomeryDivisorDataPool.Shared; + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); + cycleLength = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: true); + divisorPool.Return(divisorData); + return cycleLength; } public void GetCycleLengths(ReadOnlySpan divisors, Span cycles) @@ -185,12 +189,14 @@ public sealed class DivisorCycleCache [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeCyclesCpu(ReadOnlySpan divisors, Span cycles, ReadOnlySpan indices) { + Queue divisorPool = MontgomeryDivisorDataPool.Shared; for (int i = 0; i < indices.Length; i++) { int targetIndex = indices[i]; ulong divisor = divisors[targetIndex]; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData); + divisorPool.Return(divisorData); } } diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index f39802a..71e8a71 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs @@ -69,6 +69,7 @@ public sealed class PrimeOrderCalculatorAccelerator AcceleratorStream stream = accelerator.CreateStream(); // TODO: Review which tables are needed when final execution path is defined. LastDigitGpuTables.WarmUp(i, stream); + _kernels[i] = new CalculatorKernels(accelerator); // HeuristicCombinedGpuTables.WarmUp(i, stream); // SmallPrimeFactorGpuTables.WarmUp(i, stream); // SharedHeuristicGpuTables.EnsureStaticTables(accelerator, stream); @@ -115,18 +116,41 @@ public sealed class PrimeOrderCalculatorAccelerator public Dictionary Pow2ModEntriesToTestOnHost = []; + private static readonly CalculatorKernels[] _kernels =new CalculatorKernels[AcceleratorPool.Shared.Accelerators.Length]; + + private sealed class CalculatorKernels(Accelerator accelerator) + { + public readonly Kernel CheckFactorsKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, Stride1D.Dense>, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(PrimeOrderGpuHeuristics.CheckFactorsKernel)); + + public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> ConvertToStandardKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelConvertToStandard)).CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); + + public readonly Action, ulong, ulong, HeuristicCombinedGpuViews> HeuristicCombinedTrialDivisionKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, HeuristicCombinedGpuViews>(PrimeTesterKernels.HeuristicTrialCombinedDivisionKernel)).CreateLauncherDelegate, ulong, ulong, HeuristicCombinedGpuViews>>(); + + public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> KeepMontgomeryKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery)).CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); + + public readonly Kernel PollardRhoKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, ArrayView1D, ArrayView1D>(Pow2MontgomeryKernels.TryPollardRhoKernel)); + + public readonly Action, GpuUInt128, ArrayView1D> Pow2ModWideKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(PrimeOrderGpuHeuristics.Pow2ModKernelWide)).CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); + + public readonly Action, ArrayView> SharesFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)).CreateLauncherDelegate, ArrayView>>(); + + public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan)).CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + + public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> SmallPrimeSieveKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel)).CreateLauncherDelegate, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>(); + + public readonly Action, int, ulong, ArrayView1D, ArrayView1D> SpecialMaxKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, int, ulong, ArrayView1D, ArrayView1D>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel)).CreateLauncherDelegate, int, ulong, ArrayView1D, ArrayView1D>>(); + } + public readonly Kernel CheckFactorsKernel; - - public readonly Action, MontgomeryDivisorData, ArrayView1D> ConvertToStandardKernelLauncher; - + public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> ConvertToStandardKernelLauncher; public readonly Action, ulong, ulong, HeuristicCombinedGpuViews> HeuristicCombinedTrialDivisionKernelLauncher; - public readonly Action, MontgomeryDivisorData, ArrayView1D> KeepMontgomeryKernelLauncher; + public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> KeepMontgomeryKernelLauncher; public readonly Kernel PollardRhoKernel; public readonly Action, GpuUInt128, ArrayView1D> Pow2ModWideKernelLauncher; public readonly Action, ArrayView> SharesFactorKernelLauncher; public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernelLauncher; public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> SmallPrimeSieveKernelLauncher; - public readonly Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D> SpecialMaxKernelLauncher; + public readonly Action, int, ulong, ArrayView1D, ArrayView1D> SpecialMaxKernelLauncher; public void EnsureCapacity(int factorsCount, int primeTesterCapacity) { @@ -260,25 +284,18 @@ public sealed class PrimeOrderCalculatorAccelerator SmallPrimeFactorSquares = tables.Squares.View; } - CheckFactorsKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D>(PrimeOrderGpuHeuristics.CheckFactorsKernel)); + var kernels = _kernels[acceleratorIndex]; - ConvertToStandardKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelConvertToStandard)).CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); - - HeuristicCombinedTrialDivisionKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, HeuristicCombinedGpuViews>(PrimeTesterKernels.HeuristicTrialCombinedDivisionKernel)).CreateLauncherDelegate, ulong, ulong, HeuristicCombinedGpuViews>>(); - - KeepMontgomeryKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery)).CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); - - PollardRhoKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, ArrayView1D, ArrayView1D>(Pow2MontgomeryKernels.TryPollardRhoKernel)); - - Pow2ModWideKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(PrimeOrderGpuHeuristics.Pow2ModKernelWide)).CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); - - SharesFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)).CreateLauncherDelegate, ArrayView>>(); - - SmallPrimeFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan)).CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); - - SmallPrimeSieveKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel)).CreateLauncherDelegate, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>(); - - SpecialMaxKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel)).CreateLauncherDelegate, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>>(); + CheckFactorsKernel = kernels.CheckFactorsKernel; + ConvertToStandardKernelLauncher = kernels.ConvertToStandardKernelLauncher; + HeuristicCombinedTrialDivisionKernelLauncher = kernels.HeuristicCombinedTrialDivisionKernelLauncher; + KeepMontgomeryKernelLauncher = kernels.KeepMontgomeryKernelLauncher; + PollardRhoKernel = kernels.PollardRhoKernel; + Pow2ModWideKernelLauncher = kernels.Pow2ModWideKernelLauncher; + SharesFactorKernelLauncher = kernels.SharesFactorKernelLauncher; + SmallPrimeFactorKernelLauncher = kernels.SmallPrimeFactorKernelLauncher; + SmallPrimeSieveKernelLauncher = kernels.SmallPrimeSieveKernelLauncher; + SpecialMaxKernelLauncher = kernels.SpecialMaxKernelLauncher; } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs index 3ac3631..9e114a1 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs @@ -17,18 +17,28 @@ internal static class Pow2MontgomeryKernels public static void Pow2MontgomeryKernelKeepMontgomery( Index1D index, ArrayView1D exponents, - MontgomeryDivisorData divisor, + ulong divisorModulus, + ulong divisorNPrime, + ulong divisorMontgomeryOne, + ulong divisorMontgomeryTwo, + ulong divisorMontgomeryTwoSquared, ArrayView1D results) { + MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); results[0] = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponents[0]); } public static void Pow2MontgomeryKernelConvertToStandard( Index1D index, ArrayView1D exponents, - MontgomeryDivisorData divisor, + ulong divisorModulus, + ulong divisorNPrime, + ulong divisorMontgomeryOne, + ulong divisorMontgomeryTwo, + ulong divisorMontgomeryTwoSquared, ArrayView1D results) { + MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); results[0] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponents[0]); } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs index 9ac135f..45e69a5 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs @@ -11,12 +11,17 @@ internal static partial class PrimeOrderGpuHeuristics int count, ulong phi, ArrayView1D, Stride1D.Dense> toTest, - MontgomeryDivisorData divisorData, + ulong divisorModulus, + ulong divisorNPrime, + ulong divisorMontgomeryOne, + ulong divisorMontgomeryTwo, + ulong divisorMontgomeryTwoSquared, ArrayView1D result ) { ulong order = phi; int entryCount = count; + MontgomeryDivisorDataGpu divisorData = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); for (int i = 0; i < entryCount; i++) { @@ -155,7 +160,11 @@ internal static partial class PrimeOrderGpuHeuristics Index1D index, ulong prime, OrderKernelConfig config, - MontgomeryDivisorData divisor, + ulong divisorModulus, + ulong divisorNPrime, + ulong divisorMontgomeryOne, + ulong divisorMontgomeryTwo, + ulong divisorMontgomeryTwoSquared, ArrayView1D primes, ArrayView1D squares, int primeCount, @@ -227,6 +236,8 @@ internal static partial class PrimeOrderGpuHeuristics SortFactors(phiFactors, phiExponents, phiFactorCount); + MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); + if (TrySpecialMaxKernel(phi, prime, phiFactors, phiFactorCount, divisor)) { resultOut[0] = phi; @@ -381,7 +392,7 @@ internal static partial class PrimeOrderGpuHeuristics ulong prime, ArrayView1D factors, int factorCount, - in MontgomeryDivisorData divisor) + in MontgomeryDivisorDataGpu divisor) { for (int i = 0; i < factorCount; i++) { @@ -401,7 +412,7 @@ internal static partial class PrimeOrderGpuHeuristics return true; } - private static ulong InitializeStartingOrderKernel(ulong prime, ulong phi, in MontgomeryDivisorData divisor) + private static ulong InitializeStartingOrderKernel(ulong prime, ulong phi, in MontgomeryDivisorDataGpu divisor) { ulong order = phi; ulong residue = prime & 7UL; @@ -423,7 +434,7 @@ internal static partial class PrimeOrderGpuHeuristics ArrayView1D factors, ArrayView1D exponents, int factorCount, - in MontgomeryDivisorData divisor) + in MontgomeryDivisorDataGpu divisor) { for (int i = 0; i < factorCount; i++) { @@ -458,7 +469,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool TryConfirmOrderKernel( ulong prime, ulong order, - in MontgomeryDivisorData divisor, + in MontgomeryDivisorDataGpu divisor, uint limit, ArrayView1D primes, ArrayView1D squares, @@ -532,7 +543,7 @@ internal static partial class PrimeOrderGpuHeuristics ulong order, ulong previousOrder, byte hasPreviousOrder, - in MontgomeryDivisorData divisor, + in MontgomeryDivisorDataGpu divisor, uint limit, int maxPowChecks, ArrayView1D primes, @@ -851,7 +862,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool TryConfirmCandidateKernel( ulong prime, ulong candidate, - in MontgomeryDivisorData divisor, + in MontgomeryDivisorDataGpu divisor, uint limit, ArrayView1D primes, ArrayView1D squares, @@ -1174,7 +1185,7 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) + private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorDataGpu divisor) { return ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent) == 1UL; } @@ -1200,9 +1211,10 @@ internal static partial class PrimeOrderGpuHeuristics /// This kernel always sets the result of the corresponding element. Callers don't need to clear the output buffers. - internal static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) + internal static void Pow2ModKernel(Index1D index, ArrayView1D exponents, ulong divisorModulus, ulong divisorNPrime, ulong divisorMontgomeryOne, ulong divisorMontgomeryTwo, ulong divisorMontgomeryTwoSquared, ArrayView1D remainders) { ulong exponent = exponents[index]; + MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent); } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs index 6ea5dd0..c785b4a 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs @@ -10,7 +10,7 @@ internal static partial class PrimeOrderGpuHeuristics ulong phi, ArrayView1D factors, int factorCount, - MontgomeryDivisorData divisor, + ulong divisorModulus, ArrayView1D candidates, ArrayView1D resultOut) { @@ -57,7 +57,7 @@ internal static partial class PrimeOrderGpuHeuristics candidates[insert + 1] = current; } - var divisorPartial = new GpuDivisorPartialData(divisor.Modulus); + var divisorPartial = new GpuDivisorPartialData(divisorModulus); var stepper = new ExponentRemainderStepperGpu(divisorPartial); if (stepper.InitializeIsUnityGpu(candidates[0])) { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 2246178..0eb7cf1 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -260,7 +260,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN ArrayView1D exponentViewDevice = exponentBuffer.View; ArrayView1D hitsView = hitsBuffer.View; ArrayView1D hitIndexView = hitIndexBuffer.View; - + Queue divisorPool = MontgomeryDivisorDataPool.Shared; while (remainingCount > 0UL && !composite) { int chunkCount = chunkCountBaseline; @@ -284,7 +284,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0) { ulong divisorValue = (ulong)nextDivisor128; - MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); + MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue); ulong divisorCycle = ResolveDivisorCycle(gpu, divisorValue, prime, in montgomeryData); if (divisorCycle == prime) { @@ -294,6 +294,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN return true; } + divisorPool.Return(montgomeryData); divisorSpan[admissibleCount] = divisorValue; divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); offsetSpan[admissibleCount] = admissibleCount; diff --git a/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs b/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs new file mode 100644 index 0000000..fcb4573 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs @@ -0,0 +1,10 @@ +namespace PerfectNumbers.Core.Gpu; + +public readonly struct MontgomeryDivisorDataGpu(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong 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; +} diff --git a/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs b/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs index 31279f5..95ecff5 100644 --- a/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs +++ b/PerfectNumbers.Core/Gpu/Pow2MontgomeryGpuCalculator.cs @@ -25,7 +25,7 @@ public static class Pow2MontgomeryGpuCalculator // resultBuffer.MemSetToZero(stream); var kernelLauncher = gpu.KeepMontgomeryKernelLauncher; - kernelLauncher(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); + kernelLauncher(stream, 1, exponentBuffer.View, divisor.Modulus, divisor.NPrime, divisor.MontgomeryOne, divisor.MontgomeryTwo, divisor.MontgomeryTwoSquared, resultBuffer.View); resultBuffer.View.CopyToCPU(stream, ref result, 1); stream.Synchronize(); @@ -54,7 +54,7 @@ public static class Pow2MontgomeryGpuCalculator var kernelLauncher = gpu.ConvertToStandardKernelLauncher!; - kernelLauncher(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); + kernelLauncher(stream, 1, exponentBuffer.View, divisor.Modulus, divisor.NPrime, divisor.MontgomeryOne, divisor.MontgomeryTwo, divisor.MontgomeryTwoSquared, resultBuffer.View); resultBuffer.View.CopyToCPU(stream, ref result, 1); stream.Synchronize(); @@ -82,7 +82,7 @@ public static class Pow2MontgomeryGpuCalculator var accelerator = gpu.Accelerator; var kernelLauncher = gpu.ConvertToStandardKernelLauncher; - kernelLauncher(stream, 1, exponentBuffer.View, divisor, resultBuffer.View); + kernelLauncher(stream, 1, exponentBuffer.View, divisor.Modulus, divisor.NPrime, divisor.MontgomeryOne, divisor.MontgomeryTwo, divisor.MontgomeryTwoSquared, resultBuffer.View); resultBuffer.View.CopyToCPU(stream, ref result, 1); stream.Synchronize(); diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 014f837..ddd453d 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -398,14 +398,18 @@ internal static partial class PrimeOrderGpuHeuristics statusBuffer.View); var orderKernel = GetOrderKernel(acceleratorIndex); - var kernelLauncher = orderKernel.CreateLauncherDelegate, ArrayView1D, int, OrderKernelBuffers>>(); + var kernelLauncher = orderKernel.CreateLauncherDelegate, ArrayView1D, int, OrderKernelBuffers>>(); kernelLauncher( stream, 1, prime, kernelConfig, - divisorData, + divisorData.Modulus, + divisorData.NPrime, + divisorData.MontgomeryOne, + divisorData.MontgomeryTwo, + divisorData.MontgomeryTwoSquared, smallPrimesView, smallSquaresView, (int)smallPrimesView.Length, @@ -417,7 +421,7 @@ internal static partial class PrimeOrderGpuHeuristics resultBuffer.View.CopyToCPU(stream, ref order, 1); stream.Synchronize(); - stream.Dispose(); + AcceleratorStreamPool.Return(acceleratorIndex, stream); phiFactorBuffer.Dispose(); phiExponentBuffer.Dispose(); workFactorBuffer.Dispose(); @@ -458,7 +462,7 @@ internal static partial class PrimeOrderGpuHeuristics } var accelerator = _accelerators[acceleratorIndex]; - var loaded = accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(CalculateOrderKernel); + var loaded = accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(CalculateOrderKernel); var kernel = KernelUtil.GetKernel(loaded); pool[acceleratorIndex] = kernel; @@ -530,11 +534,12 @@ internal static partial class PrimeOrderGpuHeuristics // remainderBuffer.MemSetToZero(stream); var pow2ModKernel = GetPow2ModKernel(acceleratorIndex); - var kernelLauncher = pow2ModKernel.CreateLauncherDelegate, MontgomeryDivisorData, ArrayView1D>>(); + var kernelLauncher = pow2ModKernel.CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); try { - kernelLauncher(stream, exponents.Length, exponentBuffer.View, divisorData, remainderBuffer.View); + kernelLauncher(stream, exponents.Length, exponentBuffer.View, divisorData.Modulus, divisorData.NPrime, divisorData.MontgomeryOne, divisorData.MontgomeryTwo, divisorData.MontgomeryTwoSquared, remainderBuffer.View); + remainderBuffer.View.CopyToCPU(stream, ref MemoryMarshal.GetReference(results), exponents.Length); stream.Synchronize(); @@ -557,7 +562,7 @@ internal static partial class PrimeOrderGpuHeuristics return true; } - private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) + private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorDataGpu divisorData, Span results) { int length = exponents.Length; for (int i = 0; i < length; i++) @@ -636,7 +641,7 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorData divisorData) + private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorDataGpu divisorData) { if (modulus <= 1UL) { @@ -655,7 +660,7 @@ internal static partial class PrimeOrderGpuHeuristics } var accelerator = _accelerators[acceleratorIndex]; - var loaded = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, ArrayView1D>(Pow2ModKernel)); + var loaded = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(Pow2ModKernel)); pool[acceleratorIndex] = loaded; return loaded; diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 1f4b2e4..ef679be 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -478,7 +478,7 @@ public sealed class HeuristicCombinedPrimeTester internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) { ulong divisor = candidate.Value; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(divisor); bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); return new(candidate, divisorData, cycleLength, hasCycleHint); } diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index 0295a49..f9f8c09 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -445,7 +445,7 @@ public sealed class HeuristicPrimeTester internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) { ulong divisor = candidate.Value; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataPool.Shared.FromModulus(divisor); bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); } diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 550186c..576ae10 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -222,8 +222,10 @@ public class MersenneDivisorCycles if (divisor <= ulong.MaxValue) { ulong divisor64 = (ulong)divisor; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64); + Queue divisorPool = MontgomeryDivisorDataPool.Shared; + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor64); ulong cycle = CalculateCycleLength(divisor64, divisorData); + divisorPool.Return(divisorData); return (UInt128)cycle; } @@ -733,6 +735,7 @@ public class MersenneDivisorCycles ulong localStart = threadStart; ulong localEnd = threadEnd; + Queue divisorPool = MontgomeryDivisorDataPool.Shared; tasks[taskIndex] = Task.Run(() => { ulong rangeLength = localEnd - localStart + 1UL; @@ -758,8 +761,9 @@ public class MersenneDivisorCycles // TODO: Migrate this generation path to the shared // unrolled-hex calculator so the CPU generator matches // the GPU benchmark leader for large divisors. - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); ulong cycle = CalculateCycleLength(divisor, divisorData); + divisorPool.Return(divisorData); localCycles[localCycleIndex++] = (divisor, cycle); } diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs index 97ba2d8..8b9e5a9 100644 --- a/PerfectNumbers.Core/MontgomeryDivisorData.cs +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -1,47 +1,24 @@ -using System.Runtime.CompilerServices; - namespace PerfectNumbers.Core; -public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared) +public sealed class MontgomeryDivisorData { - public readonly ulong Modulus = modulus; + public static readonly MontgomeryDivisorData Empty = new(); + public ulong Modulus; + public ulong NPrime; + public ulong MontgomeryOne; + public ulong MontgomeryTwo; + public ulong MontgomeryTwoSquared; - public readonly ulong NPrime = nPrime; - - public readonly ulong MontgomeryOne = montgomeryOne; - - public readonly ulong MontgomeryTwo = montgomeryTwo; - public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - 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 static MontgomeryDivisorData FromModulus(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 new MontgomeryDivisorData( - modulus, - nPrime, - montgomeryOne, - montgomeryTwo, - montgomeryTwoSquared); - } + public MontgomeryDivisorData() + { + } + public MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared) + { + Modulus = modulus; + NPrime = nPrime; + MontgomeryOne = montgomeryOne; + MontgomeryTwo = montgomeryTwo; + MontgomeryTwoSquared = montgomeryTwoSquared; + } } diff --git a/PerfectNumbers.Core/MontgomeryDivisorDataPool.cs b/PerfectNumbers.Core/MontgomeryDivisorDataPool.cs new file mode 100644 index 0000000..26bf81a --- /dev/null +++ b/PerfectNumbers.Core/MontgomeryDivisorDataPool.cs @@ -0,0 +1,80 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +public static class MontgomeryDivisorDataPool +{ + [ThreadStatic] + private static Queue? _pool; + public static Queue Shared => _pool ??= new(); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + 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 static MontgomeryDivisorData Rent(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared) + { + var pool = _pool ??= new(); + if (pool.TryDequeue(out var 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)] + public static void Return(MontgomeryDivisorData divisorData) => _pool!.Enqueue(divisorData); + + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static MontgomeryDivisorData FromModulus(this Queue 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); + } + + public static MontgomeryDivisorData Rent(this Queue queue, ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared) + { + if (queue.TryDequeue(out var 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)] + public static void Return(this Queue queue, MontgomeryDivisorData divisorData) => queue.Enqueue(divisorData); +} diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index dceda8a..1f6743f 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -8,7 +8,7 @@ public static class PerfectNumberConstants public const int DefaultSmallPrimeFactorSlotCount = 64; public const int DefaultSpecialMaxFactorCapacity = 1024; public const ulong ExtraDivisorCycleSearchLimit = 64UL; - public const int GpuFrequency = 16 - 1; + public const int GpuFrequency = 4 - 1; public const int MaxQForDivisorCycles = 4_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxOddPowersCount = 128; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index a7cc4d4..7dc2454 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -1389,7 +1389,7 @@ internal static partial class PrimeOrderCalculator private static bool RunOnGpu() { // return false; - // return true; + return true; // int cpuCount; // lock(_sync) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index 3796425..1857a52 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -44,7 +44,7 @@ internal static partial class PrimeOrderCalculator // var kernelLauncher = gpu.CheckFactorsKernelLauncher; var kernel = gpu.CheckFactorsKernel; - kernel.Launch(stream, 1, entryCount, phi, gpu.Pow2ModEntriesToTestOnDevice.View, divisorData, gpu.OutputUlong.View); + kernel.Launch(stream, 1, entryCount, phi, gpu.Pow2ModEntriesToTestOnDevice.View, divisorData.Modulus, divisorData.NPrime, divisorData.MontgomeryOne, divisorData.MontgomeryTwo, divisorData.MontgomeryTwoSquared, gpu.OutputUlong.View); gpu.OutputUlong.View.CopyToCPU(stream, ref order, 1); stream.Synchronize(); @@ -146,7 +146,7 @@ internal static partial class PrimeOrderCalculator phi, specialMaxFactorsView, factorCount, - divisorData, + divisorData.Modulus, gpu.OutputUlong.View, specialMaxResultView); diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs index 2fe2791..8734c4d 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -82,7 +82,7 @@ internal static partial class PrimeOrderCalculator Console.WriteLine(messageFactory()); } } - private static UInt128 CalculateWideInternal(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + private static UInt128 CalculateWideInternal(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData? divisorData, in PrimeOrderSearchConfig config) { if (prime <= UInt128.One) { diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index d3fd7e6..70180d1 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -86,13 +86,15 @@ internal static partial class PrimeOrderCalculator } ulong prime64 = (ulong)prime; - divisorData = MontgomeryDivisorData.FromModulus(prime64); + Queue divisorPool = MontgomeryDivisorDataPool.Shared; + divisorData = divisorPool.FromModulus(prime64); ulong order64 = Calculate(prime64, previous, divisorData, config, device); + divisorPool.Return(divisorData); result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; } else { - divisorData = default; + divisorData = MontgomeryDivisorData.Empty; result = CalculateWideInternal(prime, previousOrder, divisorData, config); } diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index 14dd54b..ba698cf 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -229,13 +229,13 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorData divisor, ulong exponent) + internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorDataGpu divisor, ulong exponent) { return Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2MontgomeryModWindowedGpuConvertToStandard(in MontgomeryDivisorData divisor, ulong exponent) + internal static ulong Pow2MontgomeryModWindowedGpuConvertToStandard(in MontgomeryDivisorDataGpu divisor, ulong exponent) { ulong modulus = divisor.Modulus; ulong nPrime = divisor.NPrime; @@ -244,7 +244,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorData divisor, ulong exponent) + private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorDataGpu divisor, ulong exponent) { ulong modulus = divisor.Modulus; // Reuse the GPU-compatible MontgomeryMultiply helper highlighted in MontgomeryMultiplyBenchmarks so the accelerator path @@ -349,7 +349,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeMontgomeryOddPowerGpu(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime) + private static ulong ComputeMontgomeryOddPowerGpu(ulong exponent, in MontgomeryDivisorDataGpu divisor, ulong modulus, ulong nPrime) { ulong baseValue = divisor.MontgomeryTwo; ulong power = divisor.MontgomeryOne;