diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index 42a425c..ad983bd 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -25,7 +25,7 @@ public class PrimeTesterBenchmarks return; } - PerfectNumberConstants.RollingAccelerators = 1; + EnvironmentConfiguration.RollingAccelerators = 1; PrimeOrderCalculatorAccelerator.WarmUp(); // HeuristicGroupABPrimeTesterAccelerator.WarmUp(); diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 7df4c28..f18579d 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -18,7 +18,9 @@ internal static class Program private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; private static MersenneNumberDivisorGpuTester? _divisorTester; private static UInt128 _divisor; - private static IMersenneNumberDivisorByDivisorTester? _byDivisorTester; + private static MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder _byDivisorCpuOrderTester; + private static MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder _byDivisorHybridOrderTester; + private static MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder _byDivisorGpuOrderTester; private static Dictionary? _byDivisorPreviousResults; private static bool _limitReached; private static ulong _byDivisorStartPrime; @@ -53,8 +55,13 @@ internal static class Program Console.WriteLine("Initializing"); int gpuPrimeThreads = Math.Max(1, _cliArguments.GpuPrimeThreads); int threadCount = Math.Max(1, _cliArguments.ThreadCount); - PerfectNumberConstants.GpuRatio = _cliArguments.GpuRatio; - PerfectNumberConstants.RollingAccelerators = Math.Min(Math.Min(PerfectNumberConstants.RollingAccelerators, gpuPrimeThreads), threadCount); + + EnvironmentConfiguration.MinK = _cliArguments.MinK; + EnvironmentConfiguration.GpuBatchSize = Math.Max(1, _cliArguments.ScanBatchSize); + EnvironmentConfiguration.GpuRatio = _cliArguments.GpuRatio; + EnvironmentConfiguration.OrderDevice = _cliArguments.OrderDevice; + EnvironmentConfiguration.RollingAccelerators = Math.Min(Math.Min(PerfectNumberConstants.RollingAccelerators, gpuPrimeThreads), threadCount); + EnvironmentConfiguration.Initialize(); // NttGpuMath.GpuTransformBackend = _cliArguments.NttBackend; // NttGpuMath.ReductionMode = _cliArguments.ModReductionMode; @@ -70,7 +77,7 @@ internal static class Program PrimeTester.GpuBatchSize = gpuPrimeBatch; GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); Console.WriteLine("Warming up GPU kernels"); - PrimeTester.WarmUpGpuKernels(PerfectNumberConstants.RollingAccelerators); + PrimeTester.WarmUpGpuKernels(EnvironmentConfiguration.RollingAccelerators); // PrimeTester.WarmUpGpuKernels(gpuPrimeThreads >> 4); Console.WriteLine("Starting up threads..."); _ = UnboundedTaskScheduler.Instance; @@ -89,11 +96,10 @@ internal static class Program : (useResidue ? PrimeTransformMode.Add : PrimeTransformMode.AddPrimes); CandidatesCalculator.Configure(transformMode); ComputationDevice mersenneDevice = _cliArguments.MersenneDevice; - ComputationDevice orderDevice = _cliArguments.OrderDevice; + ComputationDevice orderDevice = EnvironmentConfiguration.OrderDevice; bool mersenneOnGpu = mersenneDevice == ComputationDevice.Gpu; - bool orderOnGpu = orderDevice == ComputationDevice.Gpu; - bool orderOnHybrid = orderDevice == ComputationDevice.Hybrid; - int scanBatchSize = Math.Max(1, _cliArguments.ScanBatchSize); + bool orderOnGpu = EnvironmentConfiguration.UseGpuOrder; + bool orderOnHybrid = EnvironmentConfiguration.UseHybridOrder; int sliceSize = Math.Max(1, _cliArguments.SliceSize); ulong residueKMax = _cliArguments.ResidueKMax; string filterFile = _cliArguments.FilterFile; @@ -235,15 +241,22 @@ internal static class Program } else if (useByDivisor) { - _byDivisorTester = mersenneOnGpu - ? new MersenneNumberDivisorByDivisorGpuTester() - : new MersenneNumberDivisorByDivisorCpuTester(orderDevice); - _byDivisorTester.MinK = _cliArguments.MinK; - _byDivisorTester.BatchSize = scanBatchSize; - if (_byDivisorTester is MersenneNumberDivisorByDivisorGpuTester && _cliArguments.MinK > ulong.MaxValue) + if (mersenneOnGpu) { - Console.WriteLine("--min-k values above UInt64 are only supported on the CPU path for --mersenne=bydivisor."); - return; + throw new NotSupportedException("--mersenne-device=gpu is not supported for --mersenne=bydivisor on this execution path."); + } + + if (EnvironmentConfiguration.UseCpuOrder) + { + _byDivisorCpuOrderTester = new MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder(); + } + else if (EnvironmentConfiguration.UseHybridOrder) + { + _byDivisorHybridOrderTester = new MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder(); + } + else + { + _byDivisorGpuOrderTester = new MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder(); } } @@ -299,7 +312,7 @@ internal static class Program NttGpuMath.GpuTransformBackend, gpuPrimeThreads, sliceSize, - scanBatchSize, + EnvironmentConfiguration.GpuBatchSize, orderWarmupLimitOverride ?? 5_000_000UL, NttGpuMath.ReductionMode, mersenneDeviceLabel, @@ -459,9 +472,11 @@ internal static class Program if (useByDivisor) { int byDivisorPrimeRange = Math.Max(1, _cliArguments.BlockSize); - MersenneNumberDivisorByDivisorTester.Run( + if (EnvironmentConfiguration.UseCpuOrder) + { + MersenneNumberDivisorByDivisorTester.Run( byDivisorCandidates, - _byDivisorTester!, + ref _byDivisorCpuOrderTester, _byDivisorPreviousResults, _byDivisorStartPrime, static () => _lastCompositeP = true, @@ -469,6 +484,33 @@ internal static class Program static (candidate, searchedMersenne, detailedCheck, passedAllTests, divisor) => CalculationResultHandler.HandleResult(candidate, divisor, searchedMersenne, detailedCheck, passedAllTests, _lastCompositeP), threadCount, byDivisorPrimeRange); + } + else if (EnvironmentConfiguration.UseHybridOrder) + { + MersenneNumberDivisorByDivisorTester.Run( + byDivisorCandidates, + ref _byDivisorHybridOrderTester, + _byDivisorPreviousResults, + _byDivisorStartPrime, + static () => _lastCompositeP = true, + static () => _lastCompositeP = false, + static (candidate, searchedMersenne, detailedCheck, passedAllTests, divisor) => CalculationResultHandler.HandleResult(candidate, divisor, searchedMersenne, detailedCheck, passedAllTests, _lastCompositeP), + threadCount, + byDivisorPrimeRange); + } + else + { + MersenneNumberDivisorByDivisorTester.Run( + byDivisorCandidates, + ref _byDivisorGpuOrderTester, + _byDivisorPreviousResults, + _byDivisorStartPrime, + static () => _lastCompositeP = true, + static () => _lastCompositeP = false, + static (candidate, searchedMersenne, detailedCheck, passedAllTests, divisor) => CalculationResultHandler.HandleResult(candidate, divisor, searchedMersenne, detailedCheck, passedAllTests, _lastCompositeP), + threadCount, + byDivisorPrimeRange); + } CalculationResultHandler.FlushBuffer(); CalculationResultHandler.ReleaseOutputBuffer(); @@ -674,9 +716,9 @@ internal static class Program // Fast residue-based composite check for p using small primes if (!_cliArguments.UseDivisor && !_cliArguments.UseByDivisor && ( - (_cliArguments.OrderDevice == ComputationDevice.Gpu && CandidatesCalculator.IsCompositeByResiduesGpu(gpu, p)) || - (_cliArguments.OrderDevice == ComputationDevice.Hybrid && CandidatesCalculator.IsCompositeByResiduesHybrid(gpu, p)) || - (_cliArguments.OrderDevice == ComputationDevice.Cpu && CandidatesCalculator.IsCompositeByResiduesCpu(p)) + (EnvironmentConfiguration.UseGpuOrder && CandidatesCalculator.IsCompositeByResiduesGpu(gpu, p)) || + (EnvironmentConfiguration.UseHybridOrder && CandidatesCalculator.IsCompositeByResiduesHybrid(gpu, p)) || + (EnvironmentConfiguration.UseCpuOrder && CandidatesCalculator.IsCompositeByResiduesCpu(p)) )) { searchedMersenne = false; @@ -711,7 +753,17 @@ internal static class Program if (_cliArguments.UseByDivisor) { // Windowed pow2mod kernel handles by-divisor scans. - return _byDivisorTester!.IsPrime(gpu, p, out detailedCheck, out _); + if (EnvironmentConfiguration.UseCpuOrder) + { + return _byDivisorCpuOrderTester.IsPrime(p, out detailedCheck, out _); + } + + if (EnvironmentConfiguration.UseHybridOrder) + { + return _byDivisorHybridOrderTester.IsPrime(p, out detailedCheck, out _); + } + + return _byDivisorGpuOrderTester.IsPrime(p, out detailedCheck, out _); } detailedCheck = MersenneTesters.Value!.IsMersennePrime(gpu, p); diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs index 9664adf..c8b4a3d 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs @@ -55,23 +55,13 @@ public class MersenneNumberDivisorByDivisorTesterTests { } - public void ResumeFromState(BigInteger lastSavedK) + public void ResumeFromState(in string stateFile, in BigInteger lastSavedK, in BigInteger minK) { } public ulong DivisorLimit => ConfiguredMaxPrime + 100UL; - public void ConfigureFromMaxPrime(ulong maxPrime) - { - ConfiguredMaxPrime = maxPrime; - } - - public ulong GetAllowedMaxDivisor(ulong prime) - { - return DivisorLimit; - } - - public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) + public bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor) { divisorsExhausted = true; if (prime % 7UL == 0UL) @@ -95,8 +85,9 @@ public class MersenneNumberDivisorByDivisorTesterTests return new DummySession(); } - public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) + public void PrepareCandidates(ulong maxPrime, in ReadOnlySpan primes, Span allowedMaxValues) { + ConfiguredMaxPrime = maxPrime; for (int i = 0; i < primes.Length; i++) { allowedMaxValues[i] = 100UL; diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index 747fa1b..0b4c3bb 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -1,125 +1,106 @@ using System.Numerics; using FluentAssertions; using PerfectNumbers.Core.Cpu; -using PerfectNumbers.Core.Gpu.Accelerators; using Xunit; namespace PerfectNumbers.Core.Tests; public class MersenneNumberDivisorCpuTesterTests { - [Fact] - [Trait("Category", "Fast")] - public void ByDivisor_tester_tracks_divisors_across_primes() - { - var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(11UL); + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_tester_tracks_divisors_across_primes() + { + var tester = new MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder(); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - try - { - tester.IsPrime(gpu, 5UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().Be(BigInteger.Zero); + ulong[] primes = [5UL, 7UL, 11UL]; + ulong[] allowedMax = new ulong[primes.Length]; + tester.PrepareCandidates(maxPrime: 11UL, primes, allowedMax); - tester.IsPrime(gpu, 7UL, out divisorsExhausted, out divisor).Should().BeTrue(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().Be(BigInteger.Zero); + tester.ResumeFromState(Path.Combine(Path.GetTempPath(), "bydivisor-cpu-tests.bin"), BigInteger.Zero, BigInteger.One); + tester.ResetStateTracking(); - tester.IsPrime(gpu, 11UL, out divisorsExhausted, out divisor).Should().BeFalse(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(BigInteger.Zero); - } - finally - { - PrimeOrderCalculatorAccelerator.Return(gpu); - } - } + tester.IsPrime(5UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); - [Fact] - [Trait("Category", "Fast")] - public void ByDivisor_session_checks_divisors_across_primes() - { - var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(13UL); + tester.IsPrime(7UL, out divisorsExhausted, out divisor).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); - try - { - ulong[] primes = [5UL, 7UL, 11UL, 13UL]; + tester.IsPrime(11UL, out divisorsExhausted, out divisor).Should().BeFalse(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().BeGreaterThan(BigInteger.Zero); + } - ulong cycle23 = MersenneDivisorCycles.CalculateCycleLengthCpu(23UL, MontgomeryDivisorData.FromModulus(23UL)); - session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes).Should().BeTrue(); + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_session_checks_divisors_across_primes() + { + var session = new MersenneCpuDivisorScanSessionWithCpuOrder(); - ulong[] primesWithoutHit = [5UL, 7UL, 13UL]; - session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primesWithoutHit).Should().BeFalse(); + ulong[] primes = [5UL, 7UL, 11UL, 13UL]; - ulong cycle31 = MersenneDivisorCycles.CalculateCycleLengthCpu(31UL, MontgomeryDivisorData.FromModulus(31UL)); - session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes).Should().BeTrue(); + ulong cycle23 = MersenneDivisorCycles.CalculateCycleLengthCpu(23UL, MontgomeryDivisorData.FromModulus(23UL)); + session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes).Should().BeTrue(); - ulong[] primesBeforeHit = [7UL, 11UL, 13UL]; - session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primesBeforeHit).Should().BeFalse(); - } - finally - { - session.Return(); - PrimeOrderCalculatorAccelerator.Return(gpu); - } - } + ulong[] primesWithoutHit = [5UL, 7UL, 13UL]; + session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primesWithoutHit).Should().BeFalse(); - [Fact] - [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite() - { - var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(43UL); + ulong cycle31 = MersenneDivisorCycles.CalculateCycleLengthCpu(31UL, MontgomeryDivisorData.FromModulus(31UL)); + session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes).Should().BeTrue(); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); - try - { - ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL]; + ulong[] primesBeforeHit = [7UL, 11UL, 13UL]; + session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primesBeforeHit).Should().BeFalse(); + } - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(7UL, divisorData); - session.CheckDivisor(7UL, divisorData, cycle, exponents).Should().BeTrue(); + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite() + { + var session = new MersenneCpuDivisorScanSessionWithCpuOrder(); - ulong[] nonDivisibleExponents = [5UL, 11UL]; - session.CheckDivisor(7UL, divisorData, cycle, nonDivisibleExponents).Should().BeFalse(); - } - finally - { - session.Return(); - PrimeOrderCalculatorAccelerator.Return(gpu); - } - } + ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL]; - [Fact] - [Trait("Category", "Fast")] - public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite() - { - var tester = new MersenneNumberDivisorByDivisorCpuTester(); - tester.ConfigureFromMaxPrime(61UL); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(7UL, divisorData); + session.CheckDivisor(7UL, divisorData, cycle, exponents).Should().BeTrue(); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); - try - { - ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL]; + ulong[] nonDivisibleExponents = [5UL, 11UL]; + session.CheckDivisor(7UL, divisorData, cycle, nonDivisibleExponents).Should().BeFalse(); + } - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); - ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(11UL, divisorData); - session.CheckDivisor(11UL, divisorData, cycle, exponents).Should().BeTrue(); + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite() + { + var session = new MersenneCpuDivisorScanSessionWithCpuOrder(); + + ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL]; + + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); + ulong cycle = MersenneDivisorCycles.CalculateCycleLengthCpu(11UL, divisorData); + session.CheckDivisor(11UL, divisorData, cycle, exponents).Should().BeTrue(); + + ulong[] nonDivisibleExponents = [5UL, 17UL]; + session.CheckDivisor(11UL, divisorData, cycle, nonDivisibleExponents).Should().BeFalse(); + } + + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_does_not_mark_mersenne_prime_as_composite_via_self_divisor() + { + var tester = new MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder(); + + tester.DivisorLimit = (BigInteger.One << 30) - BigInteger.One; + BigInteger kSelf = ((BigInteger.One << 31) - 2) / 62; // (2^p - 2) / (2p) for p=31. + + tester.ResumeFromState(Path.Combine(Path.GetTempPath(), "bydivisor-selfdivisor.bin"), BigInteger.Zero, kSelf); + tester.ResetStateTracking(); + + tester.IsPrime(31UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); + } - ulong[] nonDivisibleExponents = [5UL, 17UL]; - session.CheckDivisor(11UL, divisorData, cycle, nonDivisibleExponents).Should().BeFalse(); - } - finally - { - session.Return(); - PrimeOrderCalculatorAccelerator.Return(gpu); - } - } } - diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 3cc6a94..ab7d730 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -68,7 +68,6 @@ public class MersenneNumberDivisorGpuTesterTests { tester.IsPrime(gpu, 35UL, 31UL, 0UL, out bool divisorsExhausted).Should().BeFalse(); divisorsExhausted.Should().BeTrue(); - } finally { @@ -93,47 +92,46 @@ public class MersenneNumberDivisorGpuTesterTests } } + private static void ConfigureByDivisorTester(MersenneNumberDivisorByDivisorGpuTester tester, ulong maxPrime) + { + tester.PrepareCandidates(maxPrime, ReadOnlySpan.Empty, Span.Empty); + tester.ResumeFromState(Path.Combine(Path.GetTempPath(), "bydivisor-gpu-tests.bin"), BigInteger.Zero, BigInteger.One); + tester.ResetStateTracking(); + } + [Fact] [Trait("Category", "Fast")] public void ByDivisor_tester_tracks_divisors_across_primes() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(43UL); + ConfigureByDivisorTester(tester, maxPrime: 43UL); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - try - { - tester.IsPrime(gpu, 31UL, out bool exhausted, out BigInteger divisor).Should().BeTrue(); - exhausted.Should().BeTrue(); - divisor.Should().Be(BigInteger.Zero); + tester.IsPrime(31UL, out bool exhausted, out BigInteger divisor).Should().BeTrue(); + exhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); - tester.IsPrime(gpu, 37UL, out exhausted, out divisor).Should().BeFalse(); - exhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(BigInteger.Zero); + tester.IsPrime(37UL, out exhausted, out divisor).Should().BeFalse(); + exhausted.Should().BeTrue(); + divisor.Should().BeGreaterThan(BigInteger.Zero); - tester.IsPrime(gpu, 43UL, out exhausted, out divisor).Should().BeFalse(); - exhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(BigInteger.Zero); - } - finally - { - PrimeOrderCalculatorAccelerator.Return(gpu); - } + tester.IsPrime(43UL, out exhausted, out divisor).Should().BeFalse(); + exhausted.Should().BeTrue(); + divisor.Should().BeGreaterThan(BigInteger.Zero); } [Fact] [Trait("Category", "Fast")] - public void ByDivisor_tester_defaults_to_zero_limit_before_configuration() + public void ByDivisor_tester_defaults_to_zero_limit_before_prepare_candidates() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.DivisorLimit.Should().Be(0UL); ulong[] primes = [31UL, 37UL]; ulong[] allowed = new ulong[primes.Length]; + tester.PrepareCandidates(maxPrime: 37UL, primes, allowed); - tester.PrepareCandidates(primes, allowed); - allowed.Should().OnlyContain(value => value == 0UL); + tester.DivisorLimit.Should().BeGreaterThan(0UL); + allowed.Should().OnlyContain(value => value > 0UL); } [Fact] @@ -141,28 +139,20 @@ public class MersenneNumberDivisorGpuTesterTests public void ByDivisor_tester_respects_filter_based_limit() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(41UL); + ConfigureByDivisorTester(tester, maxPrime: 41UL); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - try - { - tester.IsPrime(gpu, 31UL, out _, out BigInteger divisor).Should().BeTrue(); - divisor.Should().Be(BigInteger.Zero); - tester.IsPrime(gpu, 37UL, out bool exhausted, out divisor).Should().BeFalse(); - exhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(BigInteger.Zero); + tester.IsPrime(31UL, out _, out BigInteger divisor).Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); - tester.IsPrime(gpu, 41UL, out bool divisorsExhausted, out divisor).Should().BeFalse(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(BigInteger.Zero); - } - finally - { - PrimeOrderCalculatorAccelerator.Return(gpu); - } + tester.IsPrime(37UL, out bool exhausted, out divisor).Should().BeFalse(); + exhausted.Should().BeTrue(); + divisor.Should().BeGreaterThan(BigInteger.Zero); + + tester.IsPrime(41UL, out bool divisorsExhausted, out divisor).Should().BeFalse(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().BeGreaterThan(BigInteger.Zero); } - [Theory] [InlineData(71UL, false)] [InlineData(89UL, true)] @@ -171,50 +161,39 @@ public class MersenneNumberDivisorGpuTesterTests public void ByDivisor_tester_matches_incremental_expectations(ulong exponent, bool expectedPrime) { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(exponent); + ConfigureByDivisorTester(tester, maxPrime: exponent); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - try + tester.IsPrime(exponent, out bool exhausted, out BigInteger divisor).Should().Be(expectedPrime); + + exhausted.Should().BeTrue(); + if (expectedPrime) { - tester.IsPrime(gpu, exponent, out bool exhausted, out BigInteger divisor).Should().Be(expectedPrime); - - exhausted.Should().BeTrue(); - if (expectedPrime) - { - divisor.Should().Be(BigInteger.Zero); - } - else - { - divisor.Should().BeGreaterThan(BigInteger.Zero); - } + divisor.Should().Be(BigInteger.Zero); } - finally + else { - PrimeOrderCalculatorAccelerator.Return(gpu); + divisor.Should().BeGreaterThan(BigInteger.Zero); } - } + [Fact] [Trait("Category", "Fast")] public void ByDivisor_gpu_tester_uses_cycle_remainders_per_divisor() { - var tester = new MersenneNumberDivisorByDivisorGpuTester - { - GpuBatchSize = 8, - }; - - tester.ConfigureFromMaxPrime(43UL); - - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); + int originalBatchSize = EnvironmentConfiguration.GpuBatchSize; try { - tester.IsPrime(gpu, 41UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeFalse(); + EnvironmentConfiguration.GpuBatchSize = 8; + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + ConfigureByDivisorTester(tester, maxPrime: 43UL); + + tester.IsPrime(41UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeFalse(); divisorsExhausted.Should().BeTrue(); divisor.Should().BeGreaterThan(BigInteger.Zero); } finally { - PrimeOrderCalculatorAccelerator.Return(gpu); + EnvironmentConfiguration.GpuBatchSize = originalBatchSize; } } @@ -223,14 +202,14 @@ public class MersenneNumberDivisorGpuTesterTests public void ByDivisor_session_computes_cycle_when_zero_is_provided() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(43UL); + ConfigureByDivisorTester(tester, maxPrime: 43UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); try { - ulong[] primes = [31UL, 37UL, 41UL]; + var session = tester.CreateDivisorSession(gpu); + ulong[] primes = [31UL, 37UL, 41UL]; session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes).Should().BeTrue(); ulong[] primesWithoutHit = [31UL, 41UL]; @@ -238,7 +217,6 @@ public class MersenneNumberDivisorGpuTesterTests } finally { - session.Return(); PrimeOrderCalculatorAccelerator.Return(gpu); } } @@ -248,12 +226,13 @@ public class MersenneNumberDivisorGpuTesterTests public void ByDivisor_session_checks_divisors_across_primes() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(43UL); + ConfigureByDivisorTester(tester, maxPrime: 43UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); try { + var session = tester.CreateDivisorSession(gpu); + ulong[] primes = [31UL, 37UL, 41UL, 43UL]; ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL)); @@ -270,7 +249,6 @@ public class MersenneNumberDivisorGpuTesterTests } finally { - session.Return(); PrimeOrderCalculatorAccelerator.Return(gpu); } } @@ -279,29 +257,34 @@ public class MersenneNumberDivisorGpuTesterTests [Trait("Category", "Fast")] public void ByDivisor_gpu_session_reuses_cycle_remainders_across_batches() { - var tester = new MersenneNumberDivisorByDivisorGpuTester - { - GpuBatchSize = 2, - }; - - tester.ConfigureFromMaxPrime(47UL); - - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); + int originalBatchSize = EnvironmentConfiguration.GpuBatchSize; try { - ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL]; + EnvironmentConfiguration.GpuBatchSize = 2; + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + ConfigureByDivisorTester(tester, maxPrime: 47UL); - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue(); + var gpu = PrimeOrderCalculatorAccelerator.Rent(1); + try + { + var session = tester.CreateDivisorSession(gpu); - ulong[] primesBeforeHit = [31UL]; - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); + ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL]; + + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue(); + + ulong[] primesBeforeHit = [31UL]; + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); + } + finally + { + PrimeOrderCalculatorAccelerator.Return(gpu); + } } finally { - session.Return(); - PrimeOrderCalculatorAccelerator.Return(gpu); + EnvironmentConfiguration.GpuBatchSize = originalBatchSize; } } @@ -309,16 +292,13 @@ public class MersenneNumberDivisorGpuTesterTests [Trait("Category", "Fast")] public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite() { - var tester = new MersenneNumberDivisorByDivisorGpuTester - { - GpuBatchSize = 5, - }; - tester.ConfigureFromMaxPrime(43UL); + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + ConfigureByDivisorTester(tester, maxPrime: 43UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); try { + var session = tester.CreateDivisorSession(gpu); ulong[] exponents = [6UL, 7UL, 9UL, 10UL, 12UL]; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL); @@ -331,7 +311,6 @@ public class MersenneNumberDivisorGpuTesterTests } finally { - session.Return(); PrimeOrderCalculatorAccelerator.Return(gpu); } } @@ -340,16 +319,13 @@ public class MersenneNumberDivisorGpuTesterTests [Trait("Category", "Fast")] public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite() { - var tester = new MersenneNumberDivisorByDivisorGpuTester - { - GpuBatchSize = 5, - }; - tester.ConfigureFromMaxPrime(61UL); + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + ConfigureByDivisorTester(tester, maxPrime: 61UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); try { + var session = tester.CreateDivisorSession(gpu); ulong[] exponents = [10UL, 11UL, 20UL, 21UL, 30UL]; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL); @@ -362,7 +338,6 @@ public class MersenneNumberDivisorGpuTesterTests } finally { - session.Return(); PrimeOrderCalculatorAccelerator.Return(gpu); } } @@ -392,13 +367,13 @@ public class MersenneNumberDivisorGpuTesterTests DivisorCycleCache.Shared.RefreshSnapshot(); var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(19UL); + ConfigureByDivisorTester(tester, maxPrime: 19UL); typeof(MersenneNumberDivisorByDivisorGpuTester) .GetField("_divisorLimit", BindingFlags.NonPublic | BindingFlags.Instance)! .SetValue(tester, 200UL); - tester.IsPrime(gpu, 19UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); + tester.IsPrime(19UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); divisorsExhausted.Should().BeTrue(); divisor.Should().Be(BigInteger.Zero); } @@ -413,32 +388,15 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_cpu_session_reuses_cycle_remainders_across_primes() + public void ByDivisor_cpu_session_checks_divisors_across_primes() { - var tester = new MersenneNumberDivisorByDivisorCpuTester - { - BatchSize = 2, - }; + var session = new MersenneCpuDivisorScanSessionWithCpuOrder(); - tester.ConfigureFromMaxPrime(47UL); + ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL]; + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthCpu(223UL, MontgomeryDivisorData.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue(); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - var session = tester.CreateDivisorSession(gpu); - try - { - ulong[] primes = [31UL, 37UL, 41UL, 43UL, 47UL]; - - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLengthGpu(gpu, 223UL, MontgomeryDivisorData.FromModulus(223UL)); - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes).Should().BeTrue(); - - ulong[] primesBeforeHit = [31UL]; - session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); - } - finally - { - session.Return(); - PrimeOrderCalculatorAccelerator.Return(gpu); - } + ulong[] primesBeforeHit = [31UL]; + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primesBeforeHit).Should().BeFalse(); } } - diff --git a/PerfectNumbers.Core/Cpu/ByteExtensions.cs b/PerfectNumbers.Core/Cpu/ByteExtensions.cs new file mode 100644 index 0000000..3720b24 --- /dev/null +++ b/PerfectNumbers.Core/Cpu/ByteExtensions.cs @@ -0,0 +1,111 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core.Cpu; + +public static class ByteExtensions +{ + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod3(this byte value, byte delta) + { + const int Modulus = 3; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod7(this byte value, byte delta) + { + const int Modulus = 7; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod13(this byte value, byte delta) + { + const int Modulus = 13; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod17(this byte value, byte delta) + { + const int Modulus = 17; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod19(this byte value, byte delta) + { + const int Modulus = 19; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod10(this byte value, byte delta) + { + const int Modulus = 10; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod11(this byte value, byte delta) + { + const int Modulus = 11; + int sum = value + delta; + + if (sum >= Modulus) + { + sum -= Modulus; + } + + return (byte)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static byte AddMod8(this byte value, byte delta) + { + return (byte)((value + delta) & 7); + } +} diff --git a/PerfectNumbers.Core/Cpu/GpuUInt128WorkSet.cs b/PerfectNumbers.Core/Cpu/GpuUInt128WorkSet.cs new file mode 100644 index 0000000..898352d --- /dev/null +++ b/PerfectNumbers.Core/Cpu/GpuUInt128WorkSet.cs @@ -0,0 +1,10 @@ +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core.Cpu; + +internal struct GpuUInt128WorkSet +{ + public GpuUInt128 Step; + public GpuUInt128 Divisor; + public GpuUInt128 Limit; +} diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs index 4dfed17..f07ca4c 100644 --- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs +++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs @@ -3,23 +3,31 @@ using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core.Cpu; -internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession -{ - private Func _getCycleLength = orderDevice switch +internal struct MersenneCpuDivisorScanSessionWithCpuOrder() : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession { - ComputationDevice.Gpu => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor, divisorData), - ComputationDevice.Hybrid => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor, divisorData), - _ => static (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthCpu(divisor, divisorData), - }; + private readonly CycleRemainderStepper _cycleStepper = new(); + private ExponentRemainderStepperCpu? _exponentStepper; - public void Configure(PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public readonly void Configure() { - _getCycleLength = orderDevice switch + } + + private static ulong _hits; + private static ulong _misses; + + private static void PrintStats(bool wasHit) + { + if (wasHit) { - ComputationDevice.Gpu => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor, divisorData), - ComputationDevice.Hybrid => (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor, divisorData), - _ => static (divisor, divisorData) => DivisorCycleCache.Shared.GetCycleLengthCpu(divisor, divisorData), - }; + Interlocked.Increment(ref _hits); + } + else + { + Interlocked.Increment(ref _misses); + } + + Console.WriteLine($"MersenneCpuDivisorScanSession, ExponentStepper hits: $_hits, misses: $_misses"); } public bool CheckDivisor( @@ -45,7 +53,154 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera if (divisorCycle == 0UL) { - divisorCycle = _getCycleLength(divisor, divisorData); + divisorCycle = DivisorCycleCache.Shared.GetCycleLengthCpu(divisor, divisorData); + + // This never kicks off in production code + // if (divisorCycle == 0UL) + // { + // // DivisorCycleCache guarantees a positive cycle for divisors greater than one. + // throw new InvalidOperationException($"Divisor cycle solver returned zero for divisor {divisor}."); + // } + } + + // 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. + CycleRemainderStepper cycleStepper = _cycleStepper; + cycleStepper.CycleLength = divisorCycle; + if (cycleStepper.Initialize(primes[0]) == 0UL) + { + return true; + } + + if (_exponentStepper is not { } exponentStepper) + { + _exponentStepper = exponentStepper = new(divisorData, divisorCycle); + } + else + { + if (!exponentStepper.MatchesDivisor(divisorData, divisorCycle)) + { + _exponentStepper = exponentStepper = new(divisorData, divisorCycle); + } + } + + if (exponentStepper.InitializeCpuIsUnity(primes[0])) + { + return true; + } + + for (int i = 1; i < length; i++) + { + if (cycleStepper.ComputeNext(primes[i]) != 0UL) + { + continue; + } + + if (exponentStepper.ComputeNextIsUnity(primes[i])) + { + return true; + } + } + + return false; + } + + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // public void Return() => ThreadStaticPools.ReturnMersenneCpuDivisorSession(this); +} + +internal struct MersenneCpuDivisorScanSessionWithHybridOrder() : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession +{ + private PrimeOrderCalculatorAccelerator? _gpu; + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void Configure(PrimeOrderCalculatorAccelerator gpu) + { + _gpu = gpu; + } + + public bool CheckDivisor( + ulong divisor, + in MontgomeryDivisorData divisorData, + ulong divisorCycle, + in ReadOnlySpan primes) + { + int length = primes.Length; + var gpu = _gpu ?? throw new InvalidOperationException("Session was not configured with an accelerator."); + if (divisorCycle == 0UL) + { + divisorCycle = DivisorCycleCache.Shared.GetCycleLengthHybrid(gpu, divisor, divisorData); + } + + // 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 cycleStepper = new CycleRemainderStepper(divisorCycle); + if (cycleStepper.Initialize(primes[0]) == 0UL) + { + return true; + } + + var exponentStepper = new ExponentRemainderStepperCpu(divisorData, divisorCycle); + if (exponentStepper.InitializeCpuIsUnity(primes[0])) + { + return true; + } + + for (int i = 1; i < length; i++) + { + if (cycleStepper.ComputeNext(primes[i]) != 0UL) + { + continue; + } + + if (exponentStepper.ComputeNextIsUnity(primes[i])) + { + return true; + } + } + + return false; + } + + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // public void Return() => ThreadStaticPools.ReturnMersenneCpuDivisorSession(this); +} + +internal struct MersenneCpuDivisorScanSessionWithGpuOrder() : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession +{ + private PrimeOrderCalculatorAccelerator? _gpu; + + public void Configure(PrimeOrderCalculatorAccelerator gpu) + { + _gpu = gpu; + } + + public bool CheckDivisor( + ulong divisor, + in MontgomeryDivisorData divisorData, + ulong divisorCycle, + in ReadOnlySpan primes) + { + int length = primes.Length; + var gpu = _gpu ?? throw new InvalidOperationException("Session was not configured with an accelerator."); + // EvenPerfectBitScanner always provides at least one exponent on this path. + // Keep the guard commented out for tests that might reuse the session differently. + // if (length == 0) + // { + // return false; + // } + + // This never kicks off in production code + // if (divisorData.Modulus != divisor) + // { + // throw new InvalidOperationException("Divisor data is for a different modulus"); + // // cachedData = MontgomeryDivisorData.FromModulus(divisor); + // } + + if (divisorCycle == 0UL) + { + divisorCycle = DivisorCycleCache.Shared.GetCycleLengthGpu(gpu, divisor, divisorData); + // This never kicks off in production code // if (divisorCycle == 0UL) // { @@ -84,6 +239,6 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera return false; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void Return() => ThreadStaticPools.ReturnMersenneCpuDivisorSession(this); + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // public void Return() => ThreadStaticPools.ReturnMersenneCpuDivisorSession(this); } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs deleted file mode 100644 index 286c3ef..0000000 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ /dev/null @@ -1,980 +0,0 @@ -using System.Numerics; -using System.Runtime.CompilerServices; -using System.Globalization; -using PerfectNumbers.Core.Gpu; -using PerfectNumbers.Core.Gpu.Accelerators; - -namespace PerfectNumbers.Core.Cpu; - -public sealed class MersenneNumberDivisorByDivisorCpuTester(ComputationDevice orderDevice = ComputationDevice.Hybrid) : IMersenneNumberDivisorByDivisorTester -{ - private BigInteger _divisorLimit; - private int _batchSize = 1_024; - private BigInteger _minK = BigInteger.One; - private readonly ComputationDevice OrderDevice = orderDevice; - private readonly TryCycleLengthDelegate _tryCalculateCycleLengthForExponent = orderDevice switch - { - ComputationDevice.Gpu => TryCalculateCycleLengthForExponentGpu, - ComputationDevice.Hybrid => TryCalculateCycleLengthForExponentHybrid, - _ => TryCalculateCycleLengthForExponentCpu, - }; - - private readonly CalculateCycleLengthDelegate _calculateCycleLength = orderDevice switch - { - ComputationDevice.Gpu => CalculateCycleLengthGpu, - ComputationDevice.Hybrid => CalculateCycleLengthHybrid, - _ => CalculateCycleLengthCpu, - }; - - private string? _stateFilePath; - private int _stateCounter; - private BigInteger _lastSavedK; - - private delegate bool TryCycleLengthDelegate( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - ulong exponent, - in MontgomeryDivisorData divisorData, - out ulong cycleLength, - out bool primeOrderFailed); - - private delegate ulong CalculateCycleLengthDelegate( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - in MontgomeryDivisorData divisorData, - bool skipPrimeOrderHeuristic); - - [ThreadStatic] - private static GpuUInt128WorkSet _divisorScanGpuWorkSet; - - private struct GpuUInt128WorkSet - { - public GpuUInt128 Step; - public GpuUInt128 Divisor; - public GpuUInt128 Limit; - } - - public int BatchSize - { - get => _batchSize; - set => _batchSize = Math.Max(1, value); - } - - public BigInteger MinK - { - get => _minK; - set => _minK = value < BigInteger.One ? BigInteger.One : value; - } - - public string? StateFilePath - { - get => _stateFilePath; - set => _stateFilePath = value; - } - - public void ResetStateTracking() - { - _stateCounter = 0; - } - - public void ResumeFromState(BigInteger lastSavedK) - { - _lastSavedK = lastSavedK; - _minK = lastSavedK + BigInteger.One; - _stateCounter = 0; - } - - public void ConfigureFromMaxPrime(ulong maxPrime) - { - _divisorLimit = ComputeDivisorLimitFromMaxPrimeBig(maxPrime); - } - - public ulong DivisorLimit - { - [MethodImpl(MethodImplOptions.AggressiveInlining)] - get - { - if (_divisorLimit.IsZero) - { - return ulong.MaxValue; - } - - return _divisorLimit > ulong.MaxValue ? ulong.MaxValue : (ulong)_divisorLimit; - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ulong GetAllowedMaxDivisor(ulong prime) => ComputeAllowedMaxDivisor(prime, _divisorLimit); - - public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) - { - BigInteger allowedMax = ComputeAllowedMaxDivisorBig(prime, _divisorLimit); - - // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. - // if (allowedMax < 3UL) - // { - // // EvenPerfectBitScanner routes primes below the small-divisor cutoff to the GPU path, so the CPU path still sees - // // trivial candidates during targeted tests. Short-circuit here to keep those runs aligned with the production flow. - // divisorsExhausted = true; - // return true; - // } - - - bool composite = CheckDivisors( - gpu, - prime, - allowedMax, - _minK, - out bool processedAll, - out divisor); - - if (composite) - { - divisorsExhausted = true; - return false; - } - - divisorsExhausted = processedAll || composite; - divisor = BigInteger.Zero; - return true; - } - - public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) - { - ulong divisorLimit = GetAllowedDivisorLimitForSpan(_divisorLimit); - int length = primes.Length; - for (int index = 0; index < length; index++) - { - allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) => ThreadStaticPools.RentMersenneCpuDivisorSession(gpu, OrderDevice); - - private bool CheckDivisors( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - BigInteger allowedMax, - BigInteger minK, - out bool processedAll, - out BigInteger foundDivisor) - { - BigInteger normalizedMinK = minK < BigInteger.One ? BigInteger.One : minK; - BigInteger step = ((BigInteger)prime) << 1; - BigInteger firstDivisor = (step * normalizedMinK) + BigInteger.One; - bool unlimited = allowedMax.IsZero; - - bool fits64 = normalizedMinK <= ulong.MaxValue && firstDivisor <= ulong.MaxValue; - bool limitFits64 = unlimited || allowedMax <= ulong.MaxValue; - if (fits64 && limitFits64) - { - ulong allowedMax64 = unlimited ? ulong.MaxValue : (ulong)allowedMax; - bool composite64 = CheckDivisors64Bit( - gpu, - prime, - allowedMax64, - (ulong)normalizedMinK, - out processedAll, - out ulong foundDivisor64); - foundDivisor = foundDivisor64; - return composite64; - } - - if (fits64 && allowedMax > ulong.MaxValue) - { - ulong allowedMax64 = ulong.MaxValue; - bool composite64 = CheckDivisors64Bit( - gpu, - prime, - allowedMax64, - (ulong)normalizedMinK, - out processedAll, - out ulong foundDivisor64); - if (composite64) - { - foundDivisor = foundDivisor64; - return true; - } - - BigInteger iterations = ((BigInteger)(allowedMax64 - (ulong)firstDivisor) / (ulong)step) + BigInteger.One; - BigInteger nextK = normalizedMinK + iterations; - return CheckDivisorsLarge( - gpu, - prime, - allowedMax, - nextK, - step, - out processedAll, - out foundDivisor); - } - - return CheckDivisorsLarge( - gpu, - prime, - allowedMax, - normalizedMinK, - step, - out processedAll, - out foundDivisor); - } - - private bool CheckDivisors64Bit( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong allowedMax, - ulong minK, - out bool processedAll, - out ulong foundDivisor) - { - foundDivisor = 0UL; - processedAll = true; - ulong currentK = minK < 1UL ? 1UL : minK; - - // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. - // Keeping the guard commented out documents the reasoning for benchmarks and tests. - // if (allowedMax < 3UL) - // { - // return false; - // } - - ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet; - - ref GpuUInt128 step = ref workSet.Step; - step.High = 0UL; - step.Low = prime; - step.ShiftLeft(1); - - ref GpuUInt128 limit = ref workSet.Limit; - limit.High = 0UL; - limit.Low = allowedMax; - - ref GpuUInt128 divisor = ref workSet.Divisor; - divisor.High = step.High; - divisor.Low = step.Low; - divisor.Add(1UL); - if (minK > 1UL) - { - GpuUInt128 offset = new(minK - 1UL); - divisor = step * offset; - divisor.Add(1UL); - } - - if (divisor.CompareTo(limit) > 0) - { - processedAll = true; - foundDivisor = 0UL; - return false; - } - - // Intentionally recomputes factorizations without a per-thread cache. - // The previous factor cache recorded virtually no hits and only slowed down the scan. - - ulong stepHigh = step.High; - ulong stepLow = step.Low; - ulong limitHigh = limit.High; - ulong divisorLow = divisor.Low; - - if (stepHigh == 0UL && limitHigh == 0UL) - { - ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL; - if (maxK == 0UL) - { - processedAll = true; - foundDivisor = 0UL; - return false; - } - - ulong startK = minK < 1UL ? 1UL : minK; - bool processedTop = true; - bool processedBottom = true; - bool composite = false; - - if (startK <= maxK) - { - composite = CheckDivisors64Range( - gpu, - prime, - stepLow, - allowedMax, - startK, - maxK, - ref currentK, - out processedTop, - out foundDivisor); - if (composite) - { - processedAll = true; - return true; - } - } - - ulong lowerEnd = startK > 1UL ? startK - 1UL : 0UL; - if (lowerEnd >= 1UL && startK == 1UL) - { - lowerEnd = Math.Min(lowerEnd, maxK); - currentK = 1UL; - composite = CheckDivisors64Range( - gpu, - prime, - stepLow, - allowedMax, - 1UL, - lowerEnd, - ref currentK, - out processedBottom, - out foundDivisor); - } - - processedAll = processedTop && processedBottom; - return composite; - } - - ulong divisorHigh = divisor.High; - var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); - - while (divisor.CompareTo(limit) <= 0) - { - if (residueStepper.IsAdmissible()) - { - ulong candidate = divisor.Low; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); - ulong divisorCycle; - // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. - if (!_tryCalculateCycleLengthForExponent( - candidate, - gpu, - 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 = _calculateCycleLength( - candidate, - gpu, - divisorData, - skipPrimeOrderHeuristic: primeOrderFailed); - } - else - { - divisorCycle = computedCycle; - } - - RecordState(currentK); - if (divisorCycle == prime) - { - // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides - // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. - foundDivisor = candidate; - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } - } - - divisor.Add(step); - currentK++; - residueStepper.Advance(); - } - - processedAll = true; - foundDivisor = 0UL; - return false; - } - - private bool CheckDivisorsLarge( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - BigInteger allowedMax, - BigInteger minK, - BigInteger step, - out bool processedAll, - out BigInteger foundDivisor) - { - processedAll = true; - foundDivisor = BigInteger.Zero; - BigInteger currentK = minK < BigInteger.One ? BigInteger.One : minK; - BigInteger divisor = (step * currentK) + BigInteger.One; - _ = gpu; - if (allowedMax.IsZero) - { - processedAll = false; - return false; - } - if (divisor > allowedMax) - { - return false; - } - - LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); - - var rem10 = new BigIntegerCycleRemainderStepper(10); - var rem8 = new BigIntegerCycleRemainderStepper(8); - var rem5 = new BigIntegerCycleRemainderStepper(5); - var rem3 = new BigIntegerCycleRemainderStepper(3); - var rem7 = new BigIntegerCycleRemainderStepper(7); - var rem11 = new BigIntegerCycleRemainderStepper(11); - var rem13 = new BigIntegerCycleRemainderStepper(13); - var rem17 = new BigIntegerCycleRemainderStepper(17); - var rem19 = new BigIntegerCycleRemainderStepper(19); - - byte remainder10 = (byte)rem10.Initialize(divisor); - byte remainder8 = (byte)rem8.Initialize(divisor); - byte remainder5 = (byte)rem5.Initialize(divisor); - byte remainder3 = (byte)rem3.Initialize(divisor); - byte remainder7 = (byte)rem7.Initialize(divisor); - byte remainder11 = (byte)rem11.Initialize(divisor); - byte remainder13 = (byte)rem13.Initialize(divisor); - byte remainder17 = (byte)rem17.Initialize(divisor); - byte remainder19 = (byte)rem19.Initialize(divisor); - - while (divisor <= allowedMax) - { - bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; - if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) - { - if (IsProbablePrimeBigInteger(divisor)) - { - BigInteger powResult = BigInteger.ModPow(2, prime, divisor); - RecordState(currentK); - if (powResult.IsOne) - { - foundDivisor = divisor; - processedAll = true; - return true; - } - } - else - { - RecordState(currentK); - } - } - - currentK += BigInteger.One; - divisor += step; - remainder10 = (byte)rem10.ComputeNext(divisor); - remainder8 = (byte)rem8.ComputeNext(divisor); - remainder5 = (byte)rem5.ComputeNext(divisor); - remainder3 = (byte)rem3.ComputeNext(divisor); - remainder7 = (byte)rem7.ComputeNext(divisor); - remainder11 = (byte)rem11.ComputeNext(divisor); - remainder13 = (byte)rem13.ComputeNext(divisor); - remainder17 = (byte)rem17.ComputeNext(divisor); - remainder19 = (byte)rem19.ComputeNext(divisor); - - if (divisor > allowedMax) - { - break; - } - } - - processedAll = true; - foundDivisor = BigInteger.Zero; - return false; - } - - private bool CheckDivisors64Range( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong step, - ulong allowedMax, - ulong startK, - ulong endK, - ref ulong currentK, - out bool processedAll, - out ulong foundDivisor) - { - if (startK < 1UL || endK < startK) - { - processedAll = true; - foundDivisor = 0UL; - return false; - } - currentK = startK; - - UInt128 step128 = step; - UInt128 allowedMax128 = allowedMax; - UInt128 startDivisor128 = (step128 * startK) + UInt128.One; - if (startDivisor128 > allowedMax128) - { - processedAll = true; - foundDivisor = 0UL; - return false; - } - - UInt128 rangeLimit128 = (step128 * endK) + UInt128.One; - if (rangeLimit128 > allowedMax128) - { - rangeLimit128 = allowedMax128; - } - - ulong startDivisor = (ulong)startDivisor128; - ulong rangeLimit = (ulong)rangeLimit128; - - var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)step128, (GpuUInt128)startDivisor128); - - return CheckDivisors64( - gpu, - prime, - step, - rangeLimit, - startDivisor, - ref residueStepper, - ref currentK, - out processedAll, - out foundDivisor); - } - -private bool CheckDivisors64( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong step, - ulong limit, - ulong divisor, - ref MersenneDivisorResidueStepper residueStepper, - ref ulong currentK, - out bool processedAll, - out ulong foundDivisor) - { - bool canAdvance = step <= limit; - if (!canAdvance) - { - if (divisor <= limit && residueStepper.IsAdmissible()) - { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - ulong divisorCycle; - if (!_tryCalculateCycleLengthForExponent( - divisor, - gpu, - prime, - divisorData, - out ulong computedCycle, - out bool primeOrderFailed) || computedCycle == 0UL) - { - divisorCycle = _calculateCycleLength( - divisor, - gpu, - divisorData, - skipPrimeOrderHeuristic: primeOrderFailed); - } - else - { - divisorCycle = computedCycle; - } - - if (divisorCycle == prime) - { - foundDivisor = divisor; - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } - } - - processedAll = true; - foundDivisor = 0UL; - return false; - } - - if (divisor > limit) - { - processedAll = true; - foundDivisor = 0UL; - return false; - } - - ulong remainingIterations = ((limit - divisor) / step) + 1UL; - while (true) - { - if (residueStepper.IsAdmissible()) - { - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - ulong divisorCycle; - if (!_tryCalculateCycleLengthForExponent( - divisor, - gpu, - prime, - divisorData, - out ulong computedCycle, - out bool primeOrderFailed) || computedCycle == 0UL) - { - divisorCycle = _calculateCycleLength( - divisor, - gpu, - divisorData, - skipPrimeOrderHeuristic: primeOrderFailed); - } - else - { - divisorCycle = computedCycle; - } - - RecordState(currentK); - if (divisorCycle == prime) - { - foundDivisor = divisor; - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } - } - - remainingIterations--; - if (remainingIterations == 0UL) - { - processedAll = true; - foundDivisor = 0UL; - return false; - } - - divisor += step; - currentK++; - residueStepper.Advance(); - } - } - - private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) - { - ulong residue = prime.Pow2MontgomeryModWithCycleConvertToStandardCpu(divisorCycle, divisorData); - return residue == 1UL ? (byte)1 : (byte)0; - } - - private void RecordState(BigInteger k) - { - string? path = _stateFilePath; - if (string.IsNullOrEmpty(path)) - { - return; - } - - if (k <= _lastSavedK) - { - return; - } - - int next = _stateCounter + 1; - if (next >= PerfectNumberConstants.ByDivisorStateSaveInterval) - { - string? directory = Path.GetDirectoryName(path); - if (!string.IsNullOrEmpty(directory)) - { - Directory.CreateDirectory(directory); - } - - File.AppendAllText(path, k.ToString(CultureInfo.InvariantCulture) + Environment.NewLine); - _stateCounter = 0; - _lastSavedK = k; - } - else - { - _stateCounter = next; - } - } - - private static bool IsProbablePrimeBigInteger(BigInteger value) - { - if (value <= 3) - { - return value >= 2; - } - - if (value.IsEven) - { - return false; - } - - BigInteger d = value - 1; - int s = 0; - while ((d & 1) == 0) - { - d >>= 1; - s++; - } - - ReadOnlySpan bases = stackalloc int[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 }; - for (int i = 0; i < bases.Length; i++) - { - int baseValue = bases[i]; - if ((BigInteger)baseValue >= value) - { - continue; - } - - var expStepper = new BigIntegerExponentStepper(value, baseValue); - BigInteger x = expStepper.Initialize(d); - if (x == 1 || x == value - 1) - { - continue; - } - - bool witnessFound = true; - for (int r = 1; r < s; r++) - { - BigInteger targetExponent = d << r; - x = expStepper.ComputeNext(targetExponent); - if (x == value - 1) - { - witnessFound = false; - break; - } - } - - if (witnessFound) - { - return false; - } - } - - return true; - } - - private static bool TryCalculateCycleLengthForExponentCpu( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - ulong exponent, - in MontgomeryDivisorData divisorData, - out ulong cycleLength, - out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( - divisor, - exponent, - divisorData, - out cycleLength, - out primeOrderFailed); - - private static bool TryCalculateCycleLengthForExponentHybrid( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - ulong exponent, - in MontgomeryDivisorData divisorData, - out ulong cycleLength, - out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentHybrid( - gpu, - divisor, - exponent, - divisorData, - out cycleLength, - out primeOrderFailed); - - private static bool TryCalculateCycleLengthForExponentGpu( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - ulong exponent, - in MontgomeryDivisorData divisorData, - out ulong cycleLength, - out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentGpu( - gpu, - divisor, - exponent, - divisorData, - out cycleLength, - out primeOrderFailed); - - private static ulong CalculateCycleLengthCpu( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - in MontgomeryDivisorData divisorData, - bool skipPrimeOrderHeuristic) => MersenneDivisorCycles.CalculateCycleLengthCpu( - divisor, - divisorData, - skipPrimeOrderHeuristic); - - private static ulong CalculateCycleLengthHybrid( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - in MontgomeryDivisorData divisorData, - bool skipPrimeOrderHeuristic) => MersenneDivisorCycles.CalculateCycleLengthHybrid( - gpu, - divisor, - divisorData, - skipPrimeOrderHeuristic); - - private static ulong CalculateCycleLengthGpu( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - in MontgomeryDivisorData divisorData, - bool skipPrimeOrderHeuristic) => MersenneDivisorCycles.CalculateCycleLengthGpu( - gpu, - divisor, - divisorData, - skipPrimeOrderHeuristic); - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod3(byte value, byte delta) - { - const int Modulus = 3; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod7(byte value, byte delta) - { - const int Modulus = 7; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod13(byte value, byte delta) - { - const int Modulus = 13; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod17(byte value, byte delta) - { - const int Modulus = 17; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod19(byte value, byte delta) - { - const int Modulus = 19; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod10(byte value, byte delta) - { - const int Modulus = 10; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod11(byte value, byte delta) - { - const int Modulus = 11; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod8(byte value, byte delta) - { - return (byte)((value + delta) & 7); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong GetAllowedDivisorLimitForSpan(BigInteger divisorLimit) - { - if (divisorLimit.IsZero) - { - return ulong.MaxValue; - } - - return divisorLimit > ulong.MaxValue ? ulong.MaxValue : (ulong)divisorLimit; - } - - private static BigInteger ComputeDivisorLimitFromMaxPrimeBig(ulong maxPrime) - { - if (maxPrime <= 1UL) - { - return BigInteger.Zero; - } - - if (maxPrime - 1UL >= 64UL) - { - return (BigInteger.One << 256) - BigInteger.One; - } - - return (BigInteger.One << (int)(maxPrime - 1UL)) - BigInteger.One; - } - - private static BigInteger ComputeAllowedMaxDivisorBig(ulong prime, BigInteger divisorLimit) - { - if (prime <= 1UL) - { - return BigInteger.Zero; - } - - if (divisorLimit.IsZero || prime - 1UL >= 64UL) - { - return divisorLimit; - } - - BigInteger candidateLimit = (BigInteger.One << (int)(prime - 1UL)) - BigInteger.One; - if (divisorLimit.IsZero) - { - return candidateLimit; - } - - return BigInteger.Min(candidateLimit, divisorLimit); - } - - private static ulong ComputeAllowedMaxDivisor(ulong prime, BigInteger divisorLimit) - { - BigInteger allowed = ComputeAllowedMaxDivisorBig(prime, divisorLimit); - if (allowed.IsZero) - { - return ulong.MaxValue; - } - - return allowed > ulong.MaxValue ? ulong.MaxValue : (ulong)allowed; - } - - -} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs new file mode 100644 index 0000000..54c0a26 --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs @@ -0,0 +1,585 @@ +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Globalization; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core.Cpu; + +public struct MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder() : IMersenneNumberDivisorByDivisorTester +{ + public BigInteger DivisorLimit; + public BigInteger MinK = EnvironmentConfiguration.MinK; + +#pragma warning disable CS8618 // StateFilePath is always set on EvenPerfectBitScanner execution path when ResumeFromState is called. + private string StateFilePath; +#pragma warning restore CS8618 // Non-nullable field must contain a non-null value when exiting constructor. Consider adding the 'required' modifier or declaring as nullable. + + private int _stateCounter; + private BigInteger _lastSavedK; + + private GpuUInt128WorkSet _divisorScanGpuWorkSet; + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void ResetStateTracking() + { + _stateCounter = 0; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void ResumeFromState(in string stateFile, in BigInteger lastSavedK, in BigInteger minK) + { + StateFilePath = stateFile; + _lastSavedK = lastSavedK; + MinK = minK; + _stateCounter = 0; + } + + public bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor) + { + BigInteger allowedMax = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisorBig(prime, DivisorLimit); + + // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. + // if (allowedMax < 3UL) + // { + // // EvenPerfectBitScanner routes primes below the small-divisor cutoff to the GPU path, so the CPU path still sees + // // trivial candidates during targeted tests. Short-circuit here to keep those runs aligned with the production flow. + // divisorsExhausted = true; + // return true; + // } + + + bool composite = CheckDivisors( + prime, + allowedMax, + MinK, + out bool processedAll, + out divisor); + + if (composite) + { + divisorsExhausted = true; + return false; + } + + divisorsExhausted = processedAll || composite; + divisor = BigInteger.Zero; + return true; + } + + public void PrepareCandidates(ulong maxPrime, in ReadOnlySpan primes, Span allowedMaxValues) + { + BigInteger divisorLimitBigInteger = MersenneNumberDivisorByDivisorTester.ComputeDivisorLimitFromMaxPrimeBig(maxPrime); + DivisorLimit = divisorLimitBigInteger; + ulong divisorLimit = MersenneNumberDivisorByDivisorTester.GetAllowedDivisorLimitForSpan(divisorLimitBigInteger); + int length = primes.Length; + for (int index = 0; index < length; index++) + { + allowedMaxValues[index] = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisor(primes[index], divisorLimit); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public readonly IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) => new MersenneCpuDivisorScanSessionWithCpuOrder(); + + private bool CheckDivisors( + ulong prime, + BigInteger allowedMax, + BigInteger minK, + out bool processedAll, + out BigInteger foundDivisor) + { + BigInteger normalizedMinK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger step = ((BigInteger)prime) << 1; + BigInteger firstDivisor = (step * normalizedMinK) + BigInteger.One; + bool unlimited = allowedMax.IsZero; + + bool fits64 = normalizedMinK <= ulong.MaxValue && firstDivisor <= ulong.MaxValue; + bool limitFits64 = unlimited || allowedMax <= ulong.MaxValue; + if (fits64 && limitFits64) + { + ulong allowedMax64 = unlimited ? ulong.MaxValue : (ulong)allowedMax; + bool composite64 = CheckDivisors64Bit( + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + foundDivisor = foundDivisor64; + return composite64; + } + + if (fits64 && allowedMax > ulong.MaxValue) + { + ulong allowedMax64 = ulong.MaxValue; + bool composite64 = CheckDivisors64Bit( + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + if (composite64) + { + foundDivisor = foundDivisor64; + return true; + } + + BigInteger iterations = ((BigInteger)(allowedMax64 - (ulong)firstDivisor) / (ulong)step) + BigInteger.One; + BigInteger nextK = normalizedMinK + iterations; + return CheckDivisorsLarge( + prime, + allowedMax, + nextK, + step, + out processedAll, + out foundDivisor); + } + + return CheckDivisorsLarge( + prime, + allowedMax, + normalizedMinK, + step, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64Bit( + ulong prime, + ulong allowedMax, + ulong minK, + out bool processedAll, + out ulong foundDivisor) + { + foundDivisor = 0UL; + ulong currentK = minK < 1UL ? 1UL : minK; + + // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. + // Keeping the guard commented out documents the reasoning for benchmarks and tests. + // if (allowedMax < 3UL) + // { + // return false; + // } + + ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet; + + ref GpuUInt128 step = ref workSet.Step; + step.High = 0UL; + step.Low = prime; + step.ShiftLeft(1); + + ref GpuUInt128 limit = ref workSet.Limit; + limit.High = 0UL; + limit.Low = allowedMax; + + ref GpuUInt128 divisor = ref workSet.Divisor; + divisor.High = step.High; + divisor.Low = step.Low; + divisor.Add(1UL); + if (minK > 1UL) + { + GpuUInt128 offset = new(minK - 1UL); + divisor = step * offset; + divisor.Add(1UL); + } + + if (divisor.CompareTo(limit) > 0) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + // Intentionally recomputes factorizations without a per-thread cache. + // The previous factor cache recorded virtually no hits and only slowed down the scan. + + ulong stepHigh = step.High; + ulong stepLow = step.Low; + ulong limitHigh = limit.High; + + if (stepHigh == 0UL && limitHigh == 0UL) + { + ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL; + if (maxK == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + ulong startK = minK < 1UL ? 1UL : minK; + bool processedTop = true; + bool processedBottom = true; + bool composite = false; + + if (startK <= maxK) + { + composite = CheckDivisors64Range( + prime, + stepLow, + allowedMax, + startK, + maxK, + ref currentK, + out processedTop, + out foundDivisor); + if (composite) + { + processedAll = true; + return true; + } + } + + ulong lowerEnd = startK > 1UL ? startK - 1UL : 0UL; + if (lowerEnd >= 1UL && startK == 1UL) + { + lowerEnd = Math.Min(lowerEnd, maxK); + currentK = 1UL; + composite = CheckDivisors64Range( + prime, + stepLow, + allowedMax, + 1UL, + lowerEnd, + ref currentK, + out processedBottom, + out foundDivisor); + } + + processedAll = processedTop && processedBottom; + return composite; + } + + var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); + + while (divisor.CompareTo(limit) <= 0) + { + if (residueStepper.IsAdmissible()) + { + ulong candidate = divisor.Low; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + ulong divisorCycle; + // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( + 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 = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( + candidate, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + RecordState(StateFilePath, currentK); + if (divisorCycle == prime) + { + // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides + // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. + foundDivisor = candidate; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + divisor.Add(step); + currentK++; + residueStepper.Advance(); + } + + processedAll = true; + foundDivisor = 0UL; + return false; + } + + private bool CheckDivisorsLarge( + ulong prime, + BigInteger allowedMax, + BigInteger minK, + BigInteger step, + out bool processedAll, + out BigInteger foundDivisor) + { + processedAll = true; + foundDivisor = BigInteger.Zero; + BigInteger currentK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger divisor = (step * currentK) + BigInteger.One; + if (allowedMax.IsZero) + { + processedAll = false; + return false; + } + if (divisor > allowedMax) + { + return false; + } + + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); + + var rem10 = new BigIntegerCycleRemainderStepper(10); + var rem8 = new BigIntegerCycleRemainderStepper(8); + var rem5 = new BigIntegerCycleRemainderStepper(5); + var rem3 = new BigIntegerCycleRemainderStepper(3); + var rem7 = new BigIntegerCycleRemainderStepper(7); + var rem11 = new BigIntegerCycleRemainderStepper(11); + var rem13 = new BigIntegerCycleRemainderStepper(13); + var rem17 = new BigIntegerCycleRemainderStepper(17); + var rem19 = new BigIntegerCycleRemainderStepper(19); + + byte remainder10 = (byte)rem10.Initialize(divisor); + byte remainder8 = (byte)rem8.Initialize(divisor); + byte remainder5 = (byte)rem5.Initialize(divisor); + byte remainder3 = (byte)rem3.Initialize(divisor); + byte remainder7 = (byte)rem7.Initialize(divisor); + byte remainder11 = (byte)rem11.Initialize(divisor); + byte remainder13 = (byte)rem13.Initialize(divisor); + byte remainder17 = (byte)rem17.Initialize(divisor); + byte remainder19 = (byte)rem19.Initialize(divisor); + + while (divisor <= allowedMax) + { + bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; + if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + { + if (MersenneNumberDivisorByDivisorTester.IsProbablePrimeBigInteger(divisor)) + { + BigInteger powResult = BigInteger.ModPow(2, prime, divisor); + RecordState(StateFilePath, currentK); + if (powResult.IsOne) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + } + else + { + RecordState(StateFilePath, currentK); + } + } + + currentK += BigInteger.One; + divisor += step; + remainder10 = (byte)rem10.ComputeNext(divisor); + remainder8 = (byte)rem8.ComputeNext(divisor); + remainder5 = (byte)rem5.ComputeNext(divisor); + remainder3 = (byte)rem3.ComputeNext(divisor); + remainder7 = (byte)rem7.ComputeNext(divisor); + remainder11 = (byte)rem11.ComputeNext(divisor); + remainder13 = (byte)rem13.ComputeNext(divisor); + remainder17 = (byte)rem17.ComputeNext(divisor); + remainder19 = (byte)rem19.ComputeNext(divisor); + + if (divisor > allowedMax) + { + break; + } + } + + processedAll = true; + foundDivisor = BigInteger.Zero; + return false; + } + + private bool CheckDivisors64Range( + ulong prime, + ulong step, + ulong allowedMax, + ulong startK, + ulong endK, + ref ulong currentK, + out bool processedAll, + out ulong foundDivisor) + { + if (startK < 1UL || endK < startK) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + currentK = startK; + + UInt128 step128 = step; + UInt128 allowedMax128 = allowedMax; + UInt128 startDivisor128 = (step128 * startK) + UInt128.One; + if (startDivisor128 > allowedMax128) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + UInt128 rangeLimit128 = (step128 * endK) + UInt128.One; + if (rangeLimit128 > allowedMax128) + { + rangeLimit128 = allowedMax128; + } + + ulong startDivisor = (ulong)startDivisor128; + ulong rangeLimit = (ulong)rangeLimit128; + + var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)step128, (GpuUInt128)startDivisor128); + + return CheckDivisors64( + prime, + step, + rangeLimit, + startDivisor, + ref residueStepper, + ref currentK, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64( + ulong prime, + ulong step, + ulong limit, + ulong divisor, + ref MersenneDivisorResidueStepper residueStepper, + ref ulong currentK, + out bool processedAll, + out ulong foundDivisor) + { + bool canAdvance = step <= limit; + if (!canAdvance) + { + if (divisor <= limit && residueStepper.IsAdmissible()) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( + divisor, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( + divisor, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + processedAll = true; + foundDivisor = 0UL; + return false; + } + + if (divisor > limit) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + ulong remainingIterations = ((limit - divisor) / step) + 1UL; + while (true) + { + if (residueStepper.IsAdmissible()) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( + divisor, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( + divisor, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + RecordState(StateFilePath, currentK); + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + remainingIterations--; + if (remainingIterations == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + divisor += step; + currentK++; + residueStepper.Advance(); + } + } + + private void RecordState(in string path, BigInteger k) + { + + if (k <= _lastSavedK) + { + return; + } + + int next = _stateCounter + 1; + if (next >= PerfectNumberConstants.ByDivisorStateSaveInterval) + { + string? directory = Path.GetDirectoryName(path); + if (!string.IsNullOrEmpty(directory)) + { + Directory.CreateDirectory(directory); + } + + File.AppendAllText(path, k.ToString(CultureInfo.InvariantCulture) + Environment.NewLine); + _stateCounter = 0; + _lastSavedK = k; + } + else + { + _stateCounter = next; + } + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder.cs new file mode 100644 index 0000000..937e3ac --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder.cs @@ -0,0 +1,600 @@ +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Globalization; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core.Cpu; + +public struct MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder() : IMersenneNumberDivisorByDivisorTester +{ + public BigInteger DivisorLimit; + public BigInteger MinK = EnvironmentConfiguration.MinK; + public readonly PrimeOrderCalculatorAccelerator Accelerator = PrimeOrderCalculatorAccelerator.Rent(1); + +#pragma warning disable CS8618 // StateFilePath is always set on EvenPerfectBitScanner execution path when ResumeFromState is called. + private string StateFilePath; +#pragma warning restore CS8618 // Non-nullable field must contain a non-null value when exiting constructor. Consider adding the 'required' modifier or declaring as nullable. + + private int _stateCounter; + private BigInteger _lastSavedK; + + private GpuUInt128WorkSet _divisorScanGpuWorkSet; + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void ResetStateTracking() + { + _stateCounter = 0; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void ResumeFromState(in string stateFile, in BigInteger lastSavedK, in BigInteger minK) + { + StateFilePath = stateFile; + _lastSavedK = lastSavedK; + MinK = minK; + _stateCounter = 0; + } + + public bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor) + { + BigInteger allowedMax = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisorBig(prime, DivisorLimit); + + // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. + // if (allowedMax < 3UL) + // { + // // EvenPerfectBitScanner routes primes below the small-divisor cutoff to the GPU path, so the CPU path still sees + // // trivial candidates during targeted tests. Short-circuit here to keep those runs aligned with the production flow. + // divisorsExhausted = true; + // return true; + // } + + + bool composite = CheckDivisors( + prime, + allowedMax, + MinK, + out bool processedAll, + out divisor); + + if (composite) + { + divisorsExhausted = true; + return false; + } + + divisorsExhausted = processedAll || composite; + divisor = BigInteger.Zero; + return true; + } + + public void PrepareCandidates(ulong maxPrime, in ReadOnlySpan primes, Span allowedMaxValues) + { + BigInteger divisorLimitBigInteger = MersenneNumberDivisorByDivisorTester.ComputeDivisorLimitFromMaxPrimeBig(maxPrime); + DivisorLimit = divisorLimitBigInteger; + ulong divisorLimit = MersenneNumberDivisorByDivisorTester.GetAllowedDivisorLimitForSpan(divisorLimitBigInteger); + int length = primes.Length; + for (int index = 0; index < length; index++) + { + allowedMaxValues[index] = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisor(primes[index], divisorLimit); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) + { + var session = new MersenneCpuDivisorScanSessionWithGpuOrder(); + session.Configure(gpu); + return session; + } + + private bool CheckDivisors( + ulong prime, + BigInteger allowedMax, + BigInteger minK, + out bool processedAll, + out BigInteger foundDivisor) + { + BigInteger normalizedMinK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger step = ((BigInteger)prime) << 1; + BigInteger firstDivisor = (step * normalizedMinK) + BigInteger.One; + bool unlimited = allowedMax.IsZero; + + bool fits64 = normalizedMinK <= ulong.MaxValue && firstDivisor <= ulong.MaxValue; + bool limitFits64 = unlimited || allowedMax <= ulong.MaxValue; + if (fits64 && limitFits64) + { + ulong allowedMax64 = unlimited ? ulong.MaxValue : (ulong)allowedMax; + bool composite64 = CheckDivisors64Bit( + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + foundDivisor = foundDivisor64; + return composite64; + } + + if (fits64 && allowedMax > ulong.MaxValue) + { + ulong allowedMax64 = ulong.MaxValue; + bool composite64 = CheckDivisors64Bit( + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + if (composite64) + { + foundDivisor = foundDivisor64; + return true; + } + + BigInteger iterations = ((BigInteger)(allowedMax64 - (ulong)firstDivisor) / (ulong)step) + BigInteger.One; + BigInteger nextK = normalizedMinK + iterations; + return CheckDivisorsLarge( + prime, + allowedMax, + nextK, + step, + out processedAll, + out foundDivisor); + } + + return CheckDivisorsLarge( + prime, + allowedMax, + normalizedMinK, + step, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64Bit( + ulong prime, + ulong allowedMax, + ulong minK, + out bool processedAll, + out ulong foundDivisor) + { + foundDivisor = 0UL; + ulong currentK = minK < 1UL ? 1UL : minK; + + // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. + // Keeping the guard commented out documents the reasoning for benchmarks and tests. + // if (allowedMax < 3UL) + // { + // return false; + // } + + ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet; + + ref GpuUInt128 step = ref workSet.Step; + step.High = 0UL; + step.Low = prime; + step.ShiftLeft(1); + + ref GpuUInt128 limit = ref workSet.Limit; + limit.High = 0UL; + limit.Low = allowedMax; + + ref GpuUInt128 divisor = ref workSet.Divisor; + divisor.High = step.High; + divisor.Low = step.Low; + divisor.Add(1UL); + if (minK > 1UL) + { + GpuUInt128 offset = new(minK - 1UL); + divisor = step * offset; + divisor.Add(1UL); + } + + if (divisor.CompareTo(limit) > 0) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + // Intentionally recomputes factorizations without a per-thread cache. + // The previous factor cache recorded virtually no hits and only slowed down the scan. + + ulong stepHigh = step.High; + ulong stepLow = step.Low; + ulong limitHigh = limit.High; + + if (stepHigh == 0UL && limitHigh == 0UL) + { + ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL; + if (maxK == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + ulong startK = minK < 1UL ? 1UL : minK; + bool processedTop = true; + bool processedBottom = true; + bool composite = false; + + if (startK <= maxK) + { + composite = CheckDivisors64Range( + prime, + stepLow, + allowedMax, + startK, + maxK, + ref currentK, + out processedTop, + out foundDivisor); + if (composite) + { + processedAll = true; + return true; + } + } + + ulong lowerEnd = startK > 1UL ? startK - 1UL : 0UL; + if (lowerEnd >= 1UL && startK == 1UL) + { + lowerEnd = Math.Min(lowerEnd, maxK); + currentK = 1UL; + composite = CheckDivisors64Range( + prime, + stepLow, + allowedMax, + 1UL, + lowerEnd, + ref currentK, + out processedBottom, + out foundDivisor); + } + + processedAll = processedTop && processedBottom; + return composite; + } + + PrimeOrderCalculatorAccelerator gpu = Accelerator; + var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); + while (divisor.CompareTo(limit) <= 0) + { + if (residueStepper.IsAdmissible()) + { + ulong candidate = divisor.Low; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + ulong divisorCycle; + // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( + candidate, + gpu, + 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 = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( + + candidate, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + RecordState(currentK); + if (divisorCycle == prime) + { + // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides + // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. + foundDivisor = candidate; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + divisor.Add(step); + currentK++; + residueStepper.Advance(); + } + + processedAll = true; + foundDivisor = 0UL; + return false; + } + + private bool CheckDivisorsLarge( + ulong prime, + BigInteger allowedMax, + BigInteger minK, + BigInteger step, + out bool processedAll, + out BigInteger foundDivisor) + { + processedAll = true; + foundDivisor = BigInteger.Zero; + BigInteger currentK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger divisor = (step * currentK) + BigInteger.One; + if (allowedMax.IsZero) + { + processedAll = false; + return false; + } + if (divisor > allowedMax) + { + return false; + } + + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); + + var rem10 = new BigIntegerCycleRemainderStepper(10); + var rem8 = new BigIntegerCycleRemainderStepper(8); + var rem5 = new BigIntegerCycleRemainderStepper(5); + var rem3 = new BigIntegerCycleRemainderStepper(3); + var rem7 = new BigIntegerCycleRemainderStepper(7); + var rem11 = new BigIntegerCycleRemainderStepper(11); + var rem13 = new BigIntegerCycleRemainderStepper(13); + var rem17 = new BigIntegerCycleRemainderStepper(17); + var rem19 = new BigIntegerCycleRemainderStepper(19); + + byte remainder10 = (byte)rem10.Initialize(divisor); + byte remainder8 = (byte)rem8.Initialize(divisor); + byte remainder5 = (byte)rem5.Initialize(divisor); + byte remainder3 = (byte)rem3.Initialize(divisor); + byte remainder7 = (byte)rem7.Initialize(divisor); + byte remainder11 = (byte)rem11.Initialize(divisor); + byte remainder13 = (byte)rem13.Initialize(divisor); + byte remainder17 = (byte)rem17.Initialize(divisor); + byte remainder19 = (byte)rem19.Initialize(divisor); + + while (divisor <= allowedMax) + { + bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; + if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + { + if (MersenneNumberDivisorByDivisorTester.IsProbablePrimeBigInteger(divisor)) + { + BigInteger powResult = BigInteger.ModPow(2, prime, divisor); + RecordState(currentK); + if (powResult.IsOne) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + } + else + { + RecordState(currentK); + } + } + + currentK += BigInteger.One; + divisor += step; + remainder10 = (byte)rem10.ComputeNext(divisor); + remainder8 = (byte)rem8.ComputeNext(divisor); + remainder5 = (byte)rem5.ComputeNext(divisor); + remainder3 = (byte)rem3.ComputeNext(divisor); + remainder7 = (byte)rem7.ComputeNext(divisor); + remainder11 = (byte)rem11.ComputeNext(divisor); + remainder13 = (byte)rem13.ComputeNext(divisor); + remainder17 = (byte)rem17.ComputeNext(divisor); + remainder19 = (byte)rem19.ComputeNext(divisor); + + if (divisor > allowedMax) + { + break; + } + } + + processedAll = true; + foundDivisor = BigInteger.Zero; + return false; + } + + private bool CheckDivisors64Range( + ulong prime, + ulong step, + ulong allowedMax, + ulong startK, + ulong endK, + ref ulong currentK, + out bool processedAll, + out ulong foundDivisor) + { + if (startK < 1UL || endK < startK) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + currentK = startK; + + UInt128 step128 = step; + UInt128 allowedMax128 = allowedMax; + UInt128 startDivisor128 = (step128 * startK) + UInt128.One; + if (startDivisor128 > allowedMax128) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + UInt128 rangeLimit128 = (step128 * endK) + UInt128.One; + if (rangeLimit128 > allowedMax128) + { + rangeLimit128 = allowedMax128; + } + + ulong startDivisor = (ulong)startDivisor128; + ulong rangeLimit = (ulong)rangeLimit128; + + var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)step128, (GpuUInt128)startDivisor128); + + return CheckDivisors64( + prime, + step, + rangeLimit, + startDivisor, + ref residueStepper, + ref currentK, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64( + ulong prime, + ulong step, + ulong limit, + ulong divisor, + ref MersenneDivisorResidueStepper residueStepper, + ref ulong currentK, + out bool processedAll, + out ulong foundDivisor) + { + PrimeOrderCalculatorAccelerator gpu = Accelerator; + + if (step > limit) + { + if (divisor <= limit && residueStepper.IsAdmissible()) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( + divisor, + gpu, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( + divisor, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + processedAll = true; + foundDivisor = 0UL; + return false; + } + + if (divisor > limit) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + ulong remainingIterations = ((limit - divisor) / step) + 1UL; + while (true) + { + if (residueStepper.IsAdmissible()) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( + divisor, + gpu, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( + divisor, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + RecordState(currentK); + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + remainingIterations--; + if (remainingIterations == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + divisor += step; + currentK++; + residueStepper.Advance(); + } + } + + private void RecordState(BigInteger k) + { + string path = StateFilePath; + + if (k <= _lastSavedK) + { + return; + } + + int next = _stateCounter + 1; + if (next >= PerfectNumberConstants.ByDivisorStateSaveInterval) + { + string? directory = Path.GetDirectoryName(path); + if (!string.IsNullOrEmpty(directory)) + { + Directory.CreateDirectory(directory); + } + + File.AppendAllText(path, k.ToString(CultureInfo.InvariantCulture) + Environment.NewLine); + _stateCounter = 0; + _lastSavedK = k; + } + else + { + _stateCounter = next; + } + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder.cs new file mode 100644 index 0000000..28f0beb --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder.cs @@ -0,0 +1,598 @@ +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Globalization; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core.Cpu; + +public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersenneNumberDivisorByDivisorTester +{ + public BigInteger DivisorLimit; + public BigInteger MinK = EnvironmentConfiguration.MinK; + public readonly PrimeOrderCalculatorAccelerator Accelerator = PrimeOrderCalculatorAccelerator.Rent(1); + +#pragma warning disable CS8618 // StateFilePath is always set on EvenPerfectBitScanner execution path when ResumeFromState is called. + private string StateFilePath; +#pragma warning restore CS8618 // Non-nullable field must contain a non-null value when exiting constructor. Consider adding the 'required' modifier or declaring as nullable. + + private int _stateCounter; + private BigInteger _lastSavedK; + + private GpuUInt128WorkSet _divisorScanGpuWorkSet; + + public void ResetStateTracking() + { + _stateCounter = 0; + } + + public void ResumeFromState(in string stateFile, in BigInteger lastSavedK, in BigInteger minK) + { + StateFilePath = stateFile; + _lastSavedK = lastSavedK; + MinK = minK; + _stateCounter = 0; + } + + public bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor) + { + BigInteger allowedMax = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisorBig(prime, DivisorLimit); + + // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. + // if (allowedMax < 3UL) + // { + // // EvenPerfectBitScanner routes primes below the small-divisor cutoff to the GPU path, so the CPU path still sees + // // trivial candidates during targeted tests. Short-circuit here to keep those runs aligned with the production flow. + // divisorsExhausted = true; + // return true; + // } + + + bool composite = CheckDivisors( + prime, + allowedMax, + MinK, + out bool processedAll, + out divisor); + + if (composite) + { + divisorsExhausted = true; + return false; + } + + divisorsExhausted = processedAll || composite; + divisor = BigInteger.Zero; + return true; + } + + public void PrepareCandidates(ulong maxPrime, in ReadOnlySpan primes, Span allowedMaxValues) + { + BigInteger divisorLimitBigInteger = MersenneNumberDivisorByDivisorTester.ComputeDivisorLimitFromMaxPrimeBig(maxPrime); + DivisorLimit = divisorLimitBigInteger; + ulong divisorLimit = MersenneNumberDivisorByDivisorTester.GetAllowedDivisorLimitForSpan(divisorLimitBigInteger); + int length = primes.Length; + for (int index = 0; index < length; index++) + { + allowedMaxValues[index] = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisor(primes[index], divisorLimit); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public readonly IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) + { + var session = new MersenneCpuDivisorScanSessionWithHybridOrder(); + session.Configure(gpu); + return session; + } + + private bool CheckDivisors( + ulong prime, + BigInteger allowedMax, + BigInteger minK, + out bool processedAll, + out BigInteger foundDivisor) + { + BigInteger normalizedMinK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger step = ((BigInteger)prime) << 1; + BigInteger firstDivisor = (step * normalizedMinK) + BigInteger.One; + bool unlimited = allowedMax.IsZero; + + bool fits64 = normalizedMinK <= ulong.MaxValue && firstDivisor <= ulong.MaxValue; + bool limitFits64 = unlimited || allowedMax <= ulong.MaxValue; + if (fits64 && limitFits64) + { + ulong allowedMax64 = unlimited ? ulong.MaxValue : (ulong)allowedMax; + bool composite64 = CheckDivisors64Bit( + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + foundDivisor = foundDivisor64; + return composite64; + } + + if (fits64 && allowedMax > ulong.MaxValue) + { + ulong allowedMax64 = ulong.MaxValue; + bool composite64 = CheckDivisors64Bit( + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + if (composite64) + { + foundDivisor = foundDivisor64; + return true; + } + + BigInteger iterations = ((BigInteger)(allowedMax64 - (ulong)firstDivisor) / (ulong)step) + BigInteger.One; + BigInteger nextK = normalizedMinK + iterations; + return CheckDivisorsLarge( + prime, + allowedMax, + nextK, + step, + out processedAll, + out foundDivisor); + } + + return CheckDivisorsLarge( + prime, + allowedMax, + normalizedMinK, + step, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64Bit( + ulong prime, + ulong allowedMax, + ulong minK, + out bool processedAll, + out ulong foundDivisor) + { + foundDivisor = 0UL; + ulong currentK = minK < 1UL ? 1UL : minK; + + // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. + // Keeping the guard commented out documents the reasoning for benchmarks and tests. + // if (allowedMax < 3UL) + // { + // return false; + // } + + ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet; + + ref GpuUInt128 step = ref workSet.Step; + step.High = 0UL; + step.Low = prime; + step.ShiftLeft(1); + + ref GpuUInt128 limit = ref workSet.Limit; + limit.High = 0UL; + limit.Low = allowedMax; + + ref GpuUInt128 divisor = ref workSet.Divisor; + divisor.High = step.High; + divisor.Low = step.Low; + divisor.Add(1UL); + if (minK > 1UL) + { + GpuUInt128 offset = new(minK - 1UL); + divisor = step * offset; + divisor.Add(1UL); + } + + if (divisor.CompareTo(limit) > 0) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + // Intentionally recomputes factorizations without a per-thread cache. + // The previous factor cache recorded virtually no hits and only slowed down the scan. + + ulong stepHigh = step.High; + ulong stepLow = step.Low; + ulong limitHigh = limit.High; + + if (stepHigh == 0UL && limitHigh == 0UL) + { + ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL; + if (maxK == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + ulong startK = minK < 1UL ? 1UL : minK; + bool processedTop = true; + bool processedBottom = true; + bool composite = false; + + if (startK <= maxK) + { + composite = CheckDivisors64Range( + prime, + stepLow, + allowedMax, + startK, + maxK, + ref currentK, + out processedTop, + out foundDivisor); + if (composite) + { + processedAll = true; + return true; + } + } + + ulong lowerEnd = startK > 1UL ? startK - 1UL : 0UL; + if (lowerEnd >= 1UL && startK == 1UL) + { + lowerEnd = Math.Min(lowerEnd, maxK); + currentK = 1UL; + composite = CheckDivisors64Range( + prime, + stepLow, + allowedMax, + 1UL, + lowerEnd, + ref currentK, + out processedBottom, + out foundDivisor); + } + + processedAll = processedTop && processedBottom; + return composite; + } + + PrimeOrderCalculatorAccelerator gpu = Accelerator; + var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); + + while (divisor.CompareTo(limit) <= 0) + { + if (residueStepper.IsAdmissible()) + { + ulong candidate = divisor.Low; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + ulong divisorCycle; + // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( + candidate, + gpu, + 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 = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( + candidate, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + RecordState(currentK); + if (divisorCycle == prime) + { + // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides + // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. + foundDivisor = candidate; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + divisor.Add(step); + currentK++; + residueStepper.Advance(); + } + + processedAll = true; + foundDivisor = 0UL; + return false; + } + + private bool CheckDivisorsLarge( + ulong prime, + BigInteger allowedMax, + BigInteger minK, + BigInteger step, + out bool processedAll, + out BigInteger foundDivisor) + { + processedAll = true; + foundDivisor = BigInteger.Zero; + BigInteger currentK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger divisor = (step * currentK) + BigInteger.One; + if (allowedMax.IsZero) + { + processedAll = false; + return false; + } + if (divisor > allowedMax) + { + return false; + } + + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); + + var rem10 = new BigIntegerCycleRemainderStepper(10); + var rem8 = new BigIntegerCycleRemainderStepper(8); + var rem5 = new BigIntegerCycleRemainderStepper(5); + var rem3 = new BigIntegerCycleRemainderStepper(3); + var rem7 = new BigIntegerCycleRemainderStepper(7); + var rem11 = new BigIntegerCycleRemainderStepper(11); + var rem13 = new BigIntegerCycleRemainderStepper(13); + var rem17 = new BigIntegerCycleRemainderStepper(17); + var rem19 = new BigIntegerCycleRemainderStepper(19); + + byte remainder10 = (byte)rem10.Initialize(divisor); + byte remainder8 = (byte)rem8.Initialize(divisor); + byte remainder5 = (byte)rem5.Initialize(divisor); + byte remainder3 = (byte)rem3.Initialize(divisor); + byte remainder7 = (byte)rem7.Initialize(divisor); + byte remainder11 = (byte)rem11.Initialize(divisor); + byte remainder13 = (byte)rem13.Initialize(divisor); + byte remainder17 = (byte)rem17.Initialize(divisor); + byte remainder19 = (byte)rem19.Initialize(divisor); + + while (divisor <= allowedMax) + { + bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; + if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + { + if (MersenneNumberDivisorByDivisorTester.IsProbablePrimeBigInteger(divisor)) + { + BigInteger powResult = BigInteger.ModPow(2, prime, divisor); + RecordState(currentK); + if (powResult.IsOne) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + } + else + { + RecordState(currentK); + } + } + + currentK += BigInteger.One; + divisor += step; + remainder10 = (byte)rem10.ComputeNext(divisor); + remainder8 = (byte)rem8.ComputeNext(divisor); + remainder5 = (byte)rem5.ComputeNext(divisor); + remainder3 = (byte)rem3.ComputeNext(divisor); + remainder7 = (byte)rem7.ComputeNext(divisor); + remainder11 = (byte)rem11.ComputeNext(divisor); + remainder13 = (byte)rem13.ComputeNext(divisor); + remainder17 = (byte)rem17.ComputeNext(divisor); + remainder19 = (byte)rem19.ComputeNext(divisor); + + if (divisor > allowedMax) + { + break; + } + } + + processedAll = true; + foundDivisor = BigInteger.Zero; + return false; + } + + private bool CheckDivisors64Range( + ulong prime, + ulong step, + ulong allowedMax, + ulong startK, + ulong endK, + ref ulong currentK, + out bool processedAll, + out ulong foundDivisor) + { + if (startK < 1UL || endK < startK) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + currentK = startK; + + UInt128 step128 = step; + UInt128 allowedMax128 = allowedMax; + UInt128 startDivisor128 = (step128 * startK) + UInt128.One; + if (startDivisor128 > allowedMax128) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + UInt128 rangeLimit128 = (step128 * endK) + UInt128.One; + if (rangeLimit128 > allowedMax128) + { + rangeLimit128 = allowedMax128; + } + + ulong startDivisor = (ulong)startDivisor128; + ulong rangeLimit = (ulong)rangeLimit128; + + var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)step128, (GpuUInt128)startDivisor128); + + return CheckDivisors64( + prime, + step, + rangeLimit, + startDivisor, + ref residueStepper, + ref currentK, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64( + ulong prime, + ulong step, + ulong limit, + ulong divisor, + ref MersenneDivisorResidueStepper residueStepper, + ref ulong currentK, + out bool processedAll, + out ulong foundDivisor) + { + PrimeOrderCalculatorAccelerator gpu = Accelerator; + bool canAdvance = step <= limit; + if (!canAdvance) + { + if (divisor <= limit && residueStepper.IsAdmissible()) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( + divisor, + gpu, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( + divisor, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + processedAll = true; + foundDivisor = 0UL; + return false; + } + + if (divisor > limit) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + ulong remainingIterations = ((limit - divisor) / step) + 1UL; + while (true) + { + if (residueStepper.IsAdmissible()) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong divisorCycle; + if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( + divisor, + gpu, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( + divisor, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + + RecordState(currentK); + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + remainingIterations--; + if (remainingIterations == 0UL) + { + processedAll = true; + foundDivisor = 0UL; + return false; + } + + divisor += step; + currentK++; + residueStepper.Advance(); + } + } + + private void RecordState(BigInteger k) + { + string path = StateFilePath; + + if (k <= _lastSavedK) + { + return; + } + + int next = _stateCounter + 1; + if (next >= PerfectNumberConstants.ByDivisorStateSaveInterval) + { + string? directory = Path.GetDirectoryName(path); + if (!string.IsNullOrEmpty(directory)) + { + Directory.CreateDirectory(directory); + } + + File.AppendAllText(path, k.ToString(CultureInfo.InvariantCulture) + Environment.NewLine); + _stateCounter = 0; + _lastSavedK = k; + } + else + { + _stateCounter = next; + } + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithCpuOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithCpuOrder.cs new file mode 100644 index 0000000..3c91308 --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithCpuOrder.cs @@ -0,0 +1,88 @@ +using System.Globalization; +using System.Numerics; +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core.Cpu; + +internal struct MersenneNumberDivisorByDivisorPrimeScanSessionWithCpuOrder( + in MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder prototype, + Action markComposite, + Action clearComposite, + Action printResult) +{ + private MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder _tester = new() + { + DivisorLimit = prototype.DivisorLimit, + MinK = EnvironmentConfiguration.MinK, + }; + + private MersenneCpuDivisorScanSessionWithCpuOrder _divisorScanSession = new(); + + private readonly Action _markComposite = markComposite; + private readonly Action _clearComposite = clearComposite; + private readonly Action _printResult = printResult; + + public void ProcessPrime(ulong prime) + { + try + { + Console.WriteLine($"Processing {prime}"); + string stateFile = Path.Combine( + PerfectNumberConstants.ByDivisorStateDirectory, + prime.ToString(CultureInfo.InvariantCulture) + ".bin"); + + ConfigureForPrime(stateFile); + + bool isPrime = _tester.IsPrime(prime, out bool divisorsExhausted, out BigInteger divisor); + + if (!isPrime) + { + if (!string.IsNullOrEmpty(stateFile) && File.Exists(stateFile)) + { + File.Delete(stateFile); + } + + _markComposite(); + _printResult(prime, true, true, false, divisor); + Console.WriteLine($"Finished processing {prime}"); + return; + } + + _clearComposite(); + _printResult(prime, true, divisorsExhausted, true, divisor); + Console.WriteLine($"Finished processing {prime}"); + } + catch (Exception ex) + { + Console.WriteLine($"Error processing {ex.Message} {ex.StackTrace}"); + Environment.Exit(1); + } + } + + private void ConfigureForPrime(string stateFile) + { + BigInteger resumeK = EnvironmentConfiguration.MinK; + if (File.Exists(stateFile)) + { + if (MersenneNumberDivisorByDivisorTester.TryReadLastSavedK(stateFile, out BigInteger lastK)) + { + resumeK = lastK + BigInteger.One; + _tester.ResumeFromState(stateFile, lastK, resumeK); + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + + _tester.ResetStateTracking(); + } + + public void Dispose() + { + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder.cs new file mode 100644 index 0000000..97a7970 --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder.cs @@ -0,0 +1,101 @@ +using System.Globalization; +using System.Numerics; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core.Cpu; + +internal struct MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder +{ + private MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder _tester; + private MersenneCpuDivisorScanSessionWithGpuOrder _divisorScanSession; + + private readonly Action _markComposite; + private readonly Action _clearComposite; + private readonly Action _printResult; + + public MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder( + in MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder prototype, + Action markComposite, + Action clearComposite, + Action printResult) + { + _tester = new MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder + { + DivisorLimit = prototype.DivisorLimit, + MinK = EnvironmentConfiguration.MinK, + }; + + _divisorScanSession = new MersenneCpuDivisorScanSessionWithGpuOrder(); + _divisorScanSession.Configure(_tester.Accelerator); + + _markComposite = markComposite; + _clearComposite = clearComposite; + _printResult = printResult; + } + + public void ProcessPrime(ulong prime) + { + try + { + Console.WriteLine($"Processing {prime}"); + string stateFile = Path.Combine( + PerfectNumberConstants.ByDivisorStateDirectory, + prime.ToString(CultureInfo.InvariantCulture) + ".bin"); + + ConfigureForPrime(stateFile); + + bool isPrime = _tester.IsPrime(prime, out bool divisorsExhausted, out BigInteger divisor); + + if (!isPrime) + { + if (!string.IsNullOrEmpty(stateFile) && File.Exists(stateFile)) + { + File.Delete(stateFile); + } + + _markComposite(); + _printResult(prime, true, true, false, divisor); + Console.WriteLine($"Finished processing {prime}"); + return; + } + + _clearComposite(); + _printResult(prime, true, divisorsExhausted, true, divisor); + Console.WriteLine($"Finished processing {prime}"); + } + catch (Exception ex) + { + Console.WriteLine($"Error processing {ex.Message} {ex.StackTrace}"); + Environment.Exit(1); + } + } + + private void ConfigureForPrime(string stateFile) + { + BigInteger resumeK = EnvironmentConfiguration.MinK; + if (File.Exists(stateFile)) + { + if (MersenneNumberDivisorByDivisorTester.TryReadLastSavedK(stateFile, out BigInteger lastK)) + { + resumeK = lastK + BigInteger.One; + _tester.ResumeFromState(stateFile, lastK, resumeK); + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + + _tester.ResetStateTracking(); + } + + public void Dispose() + { + PrimeOrderCalculatorAccelerator.Return(_tester.Accelerator); + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder.cs new file mode 100644 index 0000000..129d540 --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder.cs @@ -0,0 +1,101 @@ +using System.Globalization; +using System.Numerics; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core.Cpu; + +internal struct MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder +{ + private MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder _tester; + private MersenneCpuDivisorScanSessionWithHybridOrder _divisorScanSession; + + private readonly Action _markComposite; + private readonly Action _clearComposite; + private readonly Action _printResult; + + public MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder( + in MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder prototype, + Action markComposite, + Action clearComposite, + Action printResult) + { + _tester = new MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder + { + DivisorLimit = prototype.DivisorLimit, + MinK = EnvironmentConfiguration.MinK, + }; + + _divisorScanSession = new MersenneCpuDivisorScanSessionWithHybridOrder(); + _divisorScanSession.Configure(_tester.Accelerator); + + _markComposite = markComposite; + _clearComposite = clearComposite; + _printResult = printResult; + } + + public void ProcessPrime(ulong prime) + { + try + { + Console.WriteLine($"Processing {prime}"); + string stateFile = Path.Combine( + PerfectNumberConstants.ByDivisorStateDirectory, + prime.ToString(CultureInfo.InvariantCulture) + ".bin"); + + ConfigureForPrime(stateFile); + + bool isPrime = _tester.IsPrime(prime, out bool divisorsExhausted, out BigInteger divisor); + + if (!isPrime) + { + if (!string.IsNullOrEmpty(stateFile) && File.Exists(stateFile)) + { + File.Delete(stateFile); + } + + _markComposite(); + _printResult(prime, true, true, false, divisor); + Console.WriteLine($"Finished processing {prime}"); + return; + } + + _clearComposite(); + _printResult(prime, true, divisorsExhausted, true, divisor); + Console.WriteLine($"Finished processing {prime}"); + } + catch (Exception ex) + { + Console.WriteLine($"Error processing {ex.Message} {ex.StackTrace}"); + Environment.Exit(1); + } + } + + private void ConfigureForPrime(string stateFile) + { + BigInteger resumeK = EnvironmentConfiguration.MinK; + if (File.Exists(stateFile)) + { + if (MersenneNumberDivisorByDivisorTester.TryReadLastSavedK(stateFile, out BigInteger lastK)) + { + resumeK = lastK + BigInteger.One; + _tester.ResumeFromState(stateFile, lastK, resumeK); + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + + _tester.ResetStateTracking(); + } + + public void Dispose() + { + PrimeOrderCalculatorAccelerator.Return(_tester.Accelerator); + } +} diff --git a/PerfectNumbers.Core/CycleRemainderStepper.cs b/PerfectNumbers.Core/CycleRemainderStepper.cs index b9297c0..5c79a24 100644 --- a/PerfectNumbers.Core/CycleRemainderStepper.cs +++ b/PerfectNumbers.Core/CycleRemainderStepper.cs @@ -4,7 +4,7 @@ namespace PerfectNumbers.Core; public struct CycleRemainderStepper(ulong cycleLength) { - private readonly ulong _cycleLength = cycleLength; + public ulong CycleLength = cycleLength; private ulong _previousPrime = 0UL; private ulong _currentRemainder = 0UL; @@ -18,7 +18,7 @@ public struct CycleRemainderStepper(ulong cycleLength) [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong Initialize(ulong prime) { - ulong cycleLength = _cycleLength; + ulong cycleLength = CycleLength; ulong remainder = prime.ReduceCycleRemainder(cycleLength); _previousPrime = prime; @@ -29,7 +29,7 @@ public struct CycleRemainderStepper(ulong cycleLength) [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong ComputeNext(ulong prime) { - ulong cycleLength = _cycleLength; + ulong cycleLength = CycleLength; ulong delta = prime - _previousPrime; _previousPrime = prime; diff --git a/PerfectNumbers.Core/EnvironmentConfiguration.cs b/PerfectNumbers.Core/EnvironmentConfiguration.cs new file mode 100644 index 0000000..6738c24 --- /dev/null +++ b/PerfectNumbers.Core/EnvironmentConfiguration.cs @@ -0,0 +1,31 @@ +using System.Numerics; +using System.Runtime.CompilerServices; +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core; + +#pragma warning disable CA2211 // The fields are settable on purpose for the best performance + +public static class EnvironmentConfiguration +{ + public static ulong DivisorLimit = ulong.MaxValue; + public static int GpuBatchSize = GpuConstants.ScanBatchSize; + public static int GpuRatio = PerfectNumberConstants.GpuRatio; + public static BigInteger MinK = BigInteger.One; + public static ComputationDevice OrderDevice = ComputationDevice.Cpu; + public static int RollingAccelerators = PerfectNumberConstants.RollingAccelerators; + public static string? StateFilePath; + public static bool UseGpuOrder; + public static bool UseHybridOrder; + public static bool UseCpuOrder; + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static void Initialize() + { + UseGpuOrder = OrderDevice == ComputationDevice.Gpu; + UseHybridOrder = OrderDevice == ComputationDevice.Hybrid; + UseCpuOrder = OrderDevice == ComputationDevice.Cpu; + } +} + +#pragma warning restore CA2211 // Non-constant fields should not be visible diff --git a/PerfectNumbers.Core/Gpu/AcceleratorPool.cs b/PerfectNumbers.Core/Gpu/AcceleratorPool.cs index be07471..a7649e8 100644 --- a/PerfectNumbers.Core/Gpu/AcceleratorPool.cs +++ b/PerfectNumbers.Core/Gpu/AcceleratorPool.cs @@ -48,7 +48,7 @@ namespace PerfectNumbers.Core.Gpu } } - public static readonly AcceleratorPool Shared = new(PerfectNumberConstants.RollingAccelerators); + public static readonly AcceleratorPool Shared = new(EnvironmentConfiguration.RollingAccelerators); private readonly int Capacity; diff --git a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs index d4c9c57..63fb2f9 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs @@ -6,8 +6,8 @@ namespace PerfectNumbers.Core.Gpu.Accelerators; public static class AcceleratorStreamPool { private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; - private static readonly PollingSemaphore[] _locks = new PollingSemaphore[PerfectNumberConstants.RollingAccelerators]; - private static readonly ConcurrentFixedCapacityStack[] _streams = new ConcurrentFixedCapacityStack[PerfectNumberConstants.RollingAccelerators]; + private static readonly PollingSemaphore[] _locks = new PollingSemaphore[EnvironmentConfiguration.RollingAccelerators]; + private static readonly ConcurrentFixedCapacityStack[] _streams = new ConcurrentFixedCapacityStack[EnvironmentConfiguration.RollingAccelerators]; private static PollingSemaphore CreateLock() => new(PerfectNumberConstants.ThreadsByAccelerator, TimeSpan.FromMilliseconds(1000), 1); diff --git a/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuTables.cs b/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuTables.cs index d56f1e7..5c4f90f 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuTables.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuTables.cs @@ -8,7 +8,7 @@ internal readonly struct HeuristicCombinedGpuTables { private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; - private static readonly HeuristicCombinedGpuTables[] _sharedTables = new HeuristicCombinedGpuTables[PerfectNumberConstants.RollingAccelerators]; + private static readonly HeuristicCombinedGpuTables[] _sharedTables = new HeuristicCombinedGpuTables[EnvironmentConfiguration.RollingAccelerators]; [MethodImpl(MethodImplOptions.AggressiveOptimization)] internal static void WarmUp(int acceleratorIndex, AcceleratorStream stream) => _sharedTables[acceleratorIndex] = new(_accelerators[acceleratorIndex], stream); diff --git a/PerfectNumbers.Core/Gpu/Accelerators/HeuristicGroupABGpuTables.cs b/PerfectNumbers.Core/Gpu/Accelerators/HeuristicGroupABGpuTables.cs index 3ca4872..a33e50c 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/HeuristicGroupABGpuTables.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/HeuristicGroupABGpuTables.cs @@ -7,9 +7,9 @@ internal sealed class HeuristicGroupABGpuTables { private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; - private static readonly HeuristicGroupABGpuTables[] _sharedTables = new HeuristicGroupABGpuTables[PerfectNumberConstants.RollingAccelerators]; + private static readonly HeuristicGroupABGpuTables[] _sharedTables = new HeuristicGroupABGpuTables[EnvironmentConfiguration.RollingAccelerators]; - private static readonly SemaphoreSlim[] _locks = [..Enumerable.Range(1, PerfectNumberConstants.RollingAccelerators).Select(_ => new SemaphoreSlim(1))]; + private static readonly SemaphoreSlim[] _locks = [..Enumerable.Range(1, EnvironmentConfiguration.RollingAccelerators).Select(_ => new SemaphoreSlim(1))]; internal static HeuristicGroupABGpuTables EnsureStaticTables(int acceleratorIndex, AcceleratorStream stream) { diff --git a/PerfectNumbers.Core/Gpu/Accelerators/LastDigitGpuTables.cs b/PerfectNumbers.Core/Gpu/Accelerators/LastDigitGpuTables.cs index d9e4a3b..af3319a 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/LastDigitGpuTables.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/LastDigitGpuTables.cs @@ -7,7 +7,7 @@ namespace PerfectNumbers.Core.Gpu.Accelerators; internal readonly struct LastDigitGpuTables { private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; - private static readonly LastDigitGpuTables[] _sharedTables = new LastDigitGpuTables[PerfectNumberConstants.RollingAccelerators]; + private static readonly LastDigitGpuTables[] _sharedTables = new LastDigitGpuTables[EnvironmentConfiguration.RollingAccelerators]; [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static LastDigitGpuTables GetStaticTables(int acceleratorIndex) => _sharedTables[acceleratorIndex]; diff --git a/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs b/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs index 2b294bd..4975e1d 100644 --- a/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs +++ b/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs @@ -9,6 +9,8 @@ public static class GpuPrimeWorkLimiter public static void Acquire() => _semaphore.Wait(); public static void Release() => _semaphore.Release(); + public static int CurrentLimit => _currentLimit; + public static void SetLimit(int value) { if (value <= 0) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuPrimeScanSession.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuPrimeScanSession.cs new file mode 100644 index 0000000..27a3003 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuPrimeScanSession.cs @@ -0,0 +1,95 @@ +using System.Globalization; +using System.Numerics; + +namespace PerfectNumbers.Core.Gpu; + +internal struct MersenneNumberDivisorByDivisorGpuPrimeScanSession +{ + private readonly MersenneNumberDivisorByDivisorGpuTester _tester; + private readonly ulong _divisorLimit; + private readonly Action _markComposite; + private readonly Action _clearComposite; + private readonly Action _printResult; + + public MersenneNumberDivisorByDivisorGpuPrimeScanSession( + MersenneNumberDivisorByDivisorGpuTester prototype, + Action markComposite, + Action clearComposite, + Action printResult) + { + _tester = new MersenneNumberDivisorByDivisorGpuTester + { + MinK = EnvironmentConfiguration.MinK, + }; + _divisorLimit = prototype.DivisorLimit; + _tester.SetDivisorLimitForClone(_divisorLimit); + + _markComposite = markComposite; + _clearComposite = clearComposite; + _printResult = printResult; + } + + public void ProcessPrime(ulong prime) + { + try + { + Console.WriteLine($"Processing {prime}"); + string stateFile = Path.Combine( + PerfectNumberConstants.ByDivisorStateDirectory, + prime.ToString(CultureInfo.InvariantCulture) + ".bin"); + + ConfigureForPrime(stateFile); + + bool isPrime = _tester.IsPrime(prime, out bool divisorsExhausted, out BigInteger divisor); + + if (!isPrime) + { + if (!string.IsNullOrEmpty(stateFile) && File.Exists(stateFile)) + { + File.Delete(stateFile); + } + + _markComposite(); + _printResult(prime, true, true, false, divisor); + Console.WriteLine($"Finished processing {prime}"); + return; + } + + _clearComposite(); + _printResult(prime, true, divisorsExhausted, true, divisor); + Console.WriteLine($"Finished processing {prime}"); + } + catch (Exception ex) + { + Console.WriteLine($"Error processing {ex.Message} {ex.StackTrace}"); + Environment.Exit(1); + } + } + + private void ConfigureForPrime(string stateFile) + { + BigInteger resumeK = EnvironmentConfiguration.MinK; + if (File.Exists(stateFile)) + { + if (MersenneNumberDivisorByDivisorTester.TryReadLastSavedK(stateFile, out BigInteger lastK)) + { + resumeK = lastK + BigInteger.One; + _tester.ResumeFromState(stateFile, lastK, resumeK); + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + } + else + { + _tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + + _tester.ResetStateTracking(); + } + + public void Dispose() + { + } +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index aaa127e..499d2d1 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -14,67 +14,46 @@ namespace PerfectNumbers.Core.Gpu; /// public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDivisorByDivisorTester { - private int _gpuBatchSize = GpuConstants.ScanBatchSize; // 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. + private readonly PrimeOrderCalculatorAccelerator Accelerator = PrimeOrderCalculatorAccelerator.Rent(1); private ulong _divisorLimit; - private BigInteger _minK = BigInteger.One; + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + internal void SetDivisorLimitForClone(ulong divisorLimit) => _divisorLimit = divisorLimit; + public BigInteger MinK = BigInteger.One; + public string? StateFilePath; private readonly ConcurrentFixedCapacityStack[] _resourcePool = [.. AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentFixedCapacityStack(PerfectNumberConstants.DefaultPoolCapacity))]; private readonly ConcurrentFixedCapacityStack _sessionPool = new(PerfectNumberConstants.DefaultPoolCapacity); - public int GpuBatchSize - { - get => _gpuBatchSize; - set => _gpuBatchSize = value < 1 ? 1 : value; - } - - int IMersenneNumberDivisorByDivisorTester.BatchSize - { - get => GpuBatchSize; - set => GpuBatchSize = value; - } - - public BigInteger MinK - { - get => _minK; - set => _minK = value < BigInteger.One ? BigInteger.One : value; - } - - public string? StateFilePath { get; set; } - + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public void ResetStateTracking() { } - public void ResumeFromState(BigInteger lastSavedK) + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void ResumeFromState(in string stateFile, in BigInteger lastSavedK, in BigInteger minK) { + StateFilePath = stateFile; if (lastSavedK > ulong.MaxValue) { - throw new NotSupportedException("GPU by-divisor tester does not support MinK above UInt64."); + ThrowNotSupportedLastSavedK(); + return; // We'll never get to here } - _minK = lastSavedK + BigInteger.One; + MinK = minK; } - public void ConfigureFromMaxPrime(ulong maxPrime) - { - // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, - // so synchronization and runtime configuration guards are unnecessary here. + private static void ThrowNotSupportedLastSavedK() => throw new NotSupportedException("GPU by-divisor tester does not support MinK above UInt64."); - _divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime); - } - - public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) + public bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor) { ulong allowedMax; - int batchCapacity; // EvenPerfectBitScanner only calls into this tester after configuring it once, so we can read the cached values without locking. allowedMax = ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); - batchCapacity = _gpuBatchSize; // Production scans never shrink the divisor window below three, so this guard stays commented out. // if (allowedMax < 3UL) @@ -84,17 +63,16 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN // } bool composite; - bool coveredRange; - ulong lastProcessed; // GpuPrimeWorkLimiter(); // int acceleratorIndex = AcceleratorPool.Shared.Rent(); // var accelerator = _accelerators[acceleratorIndex]; + PrimeOrderCalculatorAccelerator gpu = Accelerator; var acceleratorIndex = gpu.AcceleratorIndex; // Monitor.Enter(gpuLease.ExecutionLock); - var resources = RentBatchResources(acceleratorIndex, batchCapacity); + var resources = RentBatchResources(acceleratorIndex, EnvironmentConfiguration.GpuBatchSize); composite = CheckDivisors( gpu, @@ -115,8 +93,8 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN resources.Offsets, resources.Counts, resources.Cycles, - out lastProcessed, - out coveredRange); + out ulong lastProcessed, + out bool coveredRange); ReturnBatchResources(acceleratorIndex, resources); // Monitor.Exit(gpuLease.ExecutionLock); @@ -134,17 +112,19 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN return true; } + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private ulong GetMinKOrThrow() { - if (_minK > ulong.MaxValue) + if (MinK > ulong.MaxValue) { - throw new NotSupportedException("GPU by-divisor tester does not support MinK above UInt64."); + ThrowNotSupportedLastSavedK(); + return 0UL; } - return (ulong)_minK; + return (ulong)MinK; } - public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) + public void PrepareCandidates(ulong maxPrime, in ReadOnlySpan primes, Span allowedMaxValues) { // The scanner always supplies matching spans, so the previous validation remains commented out. // if (allowedMaxValues.Length < primes.Length) @@ -165,7 +145,12 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN // divisorLimit = _divisorLimit; // } - divisorLimit = _divisorLimit; + // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, + // so synchronization and runtime configuration guards are unnecessary here. + + divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime); + + _divisorLimit = divisorLimit; for (int index = 0; index < primes.Length; index++) { @@ -526,12 +511,12 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN private readonly Accelerator _accelerator; private readonly PrimeOrderCalculatorAccelerator _primeOrderCalculatorAccelerator; private readonly AcceleratorStream _stream; - private Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> _kernel = null!; - private MemoryBuffer1D _divisorBuffer = null!; - private MemoryBuffer1D _offsetBuffer = null!; - private MemoryBuffer1D _countBuffer = null!; - private MemoryBuffer1D _cycleBuffer = null!; - private MemoryBuffer1D _firstHitBuffer = null!; + private readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> _kernel = null!; + private readonly MemoryBuffer1D _divisorBuffer = null!; + private readonly MemoryBuffer1D _offsetBuffer = null!; + private readonly MemoryBuffer1D _countBuffer = null!; + private readonly MemoryBuffer1D _cycleBuffer = null!; + private readonly MemoryBuffer1D _firstHitBuffer = null!; private MemoryBuffer1D _exponentsBuffer = null!; private MemoryBuffer1D _hitBuffer = null!; private ulong[] _hostBuffer = null!; @@ -546,7 +531,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN Accelerator accelerator = gpu.Accelerator; _accelerator = accelerator; _stream = accelerator.CreateStream(); - _capacity = Math.Max(1, owner._gpuBatchSize); // lock (accelerator) { @@ -559,11 +543,9 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN int desiredCapacity = PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount; _capacity = desiredCapacity; - // lock (accelerator) - { - _exponentsBuffer = accelerator.Allocate1D(desiredCapacity); - _hitBuffer = accelerator.Allocate1D(desiredCapacity); - } + + _exponentsBuffer = accelerator.Allocate1D(desiredCapacity); + _hitBuffer = accelerator.Allocate1D(desiredCapacity); _hostBuffer = ThreadStaticPools.UlongPool.Rent(desiredCapacity); @@ -729,23 +711,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); } - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - public ulong GetAllowedMaxDivisor(ulong prime) - { - // The configuration is immutable after setup, so the cached limit can be read without locking. - // lock (_sync) - // { - // if (!_isConfigured) - // { - // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - // } - - // return ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); - // } - - return ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); - } - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) { diff --git a/PerfectNumbers.Core/Gpu/SmallPrimeFactorGpuTables.cs b/PerfectNumbers.Core/Gpu/SmallPrimeFactorGpuTables.cs index 614ca98..7ceba3e 100644 --- a/PerfectNumbers.Core/Gpu/SmallPrimeFactorGpuTables.cs +++ b/PerfectNumbers.Core/Gpu/SmallPrimeFactorGpuTables.cs @@ -7,7 +7,7 @@ public readonly struct SmallPrimeFactorGpuTables { private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; - private static readonly SmallPrimeFactorGpuTables?[] _sharedTables = new SmallPrimeFactorGpuTables?[PerfectNumberConstants.RollingAccelerators]; + private static readonly SmallPrimeFactorGpuTables?[] _sharedTables = new SmallPrimeFactorGpuTables?[EnvironmentConfiguration.RollingAccelerators]; [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static SmallPrimeFactorGpuTables? GetStaticTables(int acceleratorIndex) => _sharedTables[acceleratorIndex]; diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index 03071d5..61e6b8b 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -5,32 +5,20 @@ namespace PerfectNumbers.Core; public interface IMersenneNumberDivisorByDivisorTester { - int BatchSize { get; set; } - - BigInteger MinK { get; set; } - - void ConfigureFromMaxPrime(ulong maxPrime); - - ulong DivisorLimit { get; } - - ulong GetAllowedMaxDivisor(ulong prime); - - bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor); + bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor); IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu); - string? StateFilePath { get; set; } - void ResetStateTracking(); - void ResumeFromState(BigInteger lastSavedK); + void ResumeFromState(in string stateFile, in BigInteger lastSavedK, in BigInteger minK); - void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues); + void PrepareCandidates(ulong maxPrime, in ReadOnlySpan primes, Span allowedMaxValues); public interface IDivisorScanSession { bool CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes); - void Return(); + // void Return(); } } diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index d04ce6b..a04eb06 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -5,21 +5,345 @@ using System.Runtime.InteropServices; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; using PerfectNumbers.Core.Cpu; +using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; public static class MersenneNumberDivisorByDivisorTester { public static void Run( - List candidates, - IMersenneNumberDivisorByDivisorTester tester, - Dictionary? previousResults, - ulong startPrime, - Action markComposite, - Action clearComposite, - Action printResult, - int threadCount, - int primesPerTask) + List candidates, + ref MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) + => RunStructTester( + candidates, + ref tester, + previousResults, + startPrime, + markComposite, + clearComposite, + printResult, + threadCount, + primesPerTask, + static ( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder prototype, + int workerCount, + int partitionSize, + Action mc, + Action cc, + Action pr) => + RunCpuOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + + public static void Run( + List candidates, + ref MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) + => RunStructTester( + candidates, + ref tester, + previousResults, + startPrime, + markComposite, + clearComposite, + printResult, + threadCount, + primesPerTask, + static ( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder prototype, + int workerCount, + int partitionSize, + Action mc, + Action cc, + Action pr) => + RunHybridOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + + public static void Run( + List candidates, + ref MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) + => RunStructTester( + candidates, + ref tester, + previousResults, + startPrime, + markComposite, + clearComposite, + printResult, + threadCount, + primesPerTask, + static ( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder prototype, + int workerCount, + int partitionSize, + Action mc, + Action cc, + Action pr) => + RunGpuOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + + internal static bool TryReadLastSavedK(string stateFile, out BigInteger lastSavedK) + { + lastSavedK = BigInteger.Zero; + try + { + bool any = false; + foreach (string line in File.ReadLines(stateFile)) + { + if (string.IsNullOrWhiteSpace(line)) + { + continue; + } + + if (BigInteger.TryParse(line, NumberStyles.None, CultureInfo.InvariantCulture, out BigInteger parsed)) + { + if (!any || parsed > lastSavedK) + { + lastSavedK = parsed; + } + + any = true; + } + } + + return any; + } + catch (IOException) + { + return false; + } + catch (UnauthorizedAccessException) + { + return false; + } + } + + public static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) + { + ulong residue = prime.Pow2MontgomeryModWithCycleConvertToStandardCpu(divisorCycle, divisorData); + return residue == 1UL ? (byte)1 : (byte)0; + } + + public static bool IsProbablePrimeBigInteger(BigInteger value) + { + if (value <= 3) + { + return value >= 2; + } + + if (value.IsEven) + { + return false; + } + + BigInteger d = value - 1; + int s = 0; + while ((d & 1) == 0) + { + d >>= 1; + s++; + } + + ReadOnlySpan bases = stackalloc int[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 }; + for (int i = 0; i < bases.Length; i++) + { + int baseValue = bases[i]; + if ((BigInteger)baseValue >= value) + { + continue; + } + + var expStepper = new BigIntegerExponentStepper(value, baseValue); + BigInteger x = expStepper.Initialize(d); + if (x == 1 || x == value - 1) + { + continue; + } + + bool witnessFound = true; + for (int r = 1; r < s; r++) + { + BigInteger targetExponent = d << r; + x = expStepper.ComputeNext(targetExponent); + if (x == value - 1) + { + witnessFound = false; + break; + } + } + + if (witnessFound) + { + return false; + } + } + + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong GetAllowedDivisorLimitForSpan(BigInteger divisorLimit) + { + if (divisorLimit.IsZero) + { + return ulong.MaxValue; + } + + return divisorLimit > ulong.MaxValue ? ulong.MaxValue : (ulong)divisorLimit; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static BigInteger ComputeDivisorLimitFromMaxPrimeBig(ulong maxPrime) + { + if (maxPrime <= 1UL) + { + return BigInteger.Zero; + } + + if (maxPrime - 1UL >= 64UL) + { + return (BigInteger.One << 256) - BigInteger.One; + } + + return (BigInteger.One << (int)(maxPrime - 1UL)) - BigInteger.One; + } + + public static BigInteger ComputeAllowedMaxDivisorBig(ulong prime, BigInteger divisorLimit) + { + if (prime <= 1UL) + { + return BigInteger.Zero; + } + + if (divisorLimit.IsZero || prime - 1UL >= 64UL) + { + return divisorLimit; + } + + BigInteger candidateLimit = (BigInteger.One << (int)(prime - 1UL)) - BigInteger.One; + if (divisorLimit.IsZero) + { + return candidateLimit; + } + + return BigInteger.Min(candidateLimit, divisorLimit); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static ulong ComputeAllowedMaxDivisor(ulong prime, BigInteger divisorLimit) + { + BigInteger allowed = ComputeAllowedMaxDivisorBig(prime, divisorLimit); + if (allowed.IsZero) + { + return ulong.MaxValue; + } + + return allowed > ulong.MaxValue ? ulong.MaxValue : (ulong)allowed; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static bool TryCalculateCycleLengthForExponentCpu( + ulong divisor, + ulong exponent, + in MontgomeryDivisorData divisorData, + out ulong cycleLength, + out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + divisor, + exponent, + divisorData, + out cycleLength, + out primeOrderFailed); + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static bool TryCalculateCycleLengthForExponentHybrid( + ulong divisor, + PrimeOrderCalculatorAccelerator gpu, + ulong exponent, + in MontgomeryDivisorData divisorData, + out ulong cycleLength, + out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentHybrid( + gpu, + divisor, + exponent, + divisorData, + out cycleLength, + out primeOrderFailed); + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static bool TryCalculateCycleLengthForExponentGpu( + ulong divisor, + PrimeOrderCalculatorAccelerator gpu, + ulong exponent, + in MontgomeryDivisorData divisorData, + out ulong cycleLength, + out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentGpu( + gpu, + divisor, + exponent, + divisorData, + out cycleLength, + out primeOrderFailed); + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static ulong CalculateCycleLengthCpu( + ulong divisor, + in MontgomeryDivisorData divisorData, + bool skipPrimeOrderHeuristic) => MersenneDivisorCycles.CalculateCycleLengthCpu( + divisor, + divisorData, + skipPrimeOrderHeuristic); + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static ulong CalculateCycleLengthHybrid( + ulong divisor, + PrimeOrderCalculatorAccelerator gpu, + in MontgomeryDivisorData divisorData, + bool skipPrimeOrderHeuristic) => MersenneDivisorCycles.CalculateCycleLengthHybrid( + gpu, + divisor, + divisorData, + skipPrimeOrderHeuristic); + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static ulong CalculateCycleLengthGpu( + ulong divisor, + PrimeOrderCalculatorAccelerator gpu, + in MontgomeryDivisorData divisorData, + bool skipPrimeOrderHeuristic) => MersenneDivisorCycles.CalculateCycleLengthGpu( + gpu, + divisor, + divisorData, + skipPrimeOrderHeuristic); + + public static void Run( + List candidates, + IMersenneNumberDivisorByDivisorTester tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) { if (candidates.Count == 0) { @@ -185,22 +509,41 @@ public static class MersenneNumberDivisorByDivisorTester // } // The by-divisor CPU configuration only feeds primes well above 1, so this fallback never executes in production runs. - tester.ConfigureFromMaxPrime(maxPrime); - int primeCount = primesToTest.Count; ulong[] primeBatch = pool.Rent(primeCount); ulong[] allowedMaxBatch = pool.Rent(primeCount); List filteredPrimes = new(primeCount); Span primeSpan = primeBatch.AsSpan(0, primeCount); + primesToTest.CopyTo(primeSpan); Span allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount); - for (int i = 0; i < primeCount; i++) - { - primeSpan[i] = primesToTest[i]; - } + // for (int i = 0; i < primeCount; i++) + // { + // primeSpan[i] = primesToTest[i]; + // } - tester.PrepareCandidates(primeSpan, allowedMaxSpan); + switch (tester) + { + case MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder typedCpu: + typedCpu.PrepareCandidates(maxPrime, primeSpan, allowedMaxSpan); + tester = typedCpu; + break; + + case MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder typedHybrid: + typedHybrid.PrepareCandidates(maxPrime, primeSpan, allowedMaxSpan); + tester = typedHybrid; + break; + + case MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder typedGpuOrder: + typedGpuOrder.PrepareCandidates(maxPrime, primeSpan, allowedMaxSpan); + tester = typedGpuOrder; + break; + + default: + tester.PrepareCandidates(maxPrime, primeSpan, allowedMaxSpan); + break; + } for (int i = 0; i < primeCount; i++) { @@ -233,7 +576,7 @@ public static class MersenneNumberDivisorByDivisorTester // return; // } - int workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount; + int workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount; // Environment.ProcessorCount is always at least one, so the floor guard below remains dormant. // if (workerCount < 1) // { @@ -241,62 +584,67 @@ public static class MersenneNumberDivisorByDivisorTester // } int chunkSize = primesPerTask <= 0 ? 1 : primesPerTask; + int partitionSize = chunkSize < 1 ? 1 : chunkSize; - IMersenneNumberDivisorByDivisorTester CreateTester(IMersenneNumberDivisorByDivisorTester prototype) + if (tester is MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder cpuOrderPrototype) { - IMersenneNumberDivisorByDivisorTester tester = prototype switch - { - MersenneNumberDivisorByDivisorCpuTester cpu => new MersenneNumberDivisorByDivisorCpuTester(cpu.OrderDevice) - { - BatchSize = prototype.BatchSize, - MinK = prototype.MinK, - }, - _ => CreateGpuTester(prototype), - }; - - tester.ConfigureFromMaxPrime(maxPrime); - return tester; + RunCpuOrderFiltered(filteredPrimes, cpuOrderPrototype, workerCount, partitionSize, markComposite, clearComposite, printResult); + return; } - static IMersenneNumberDivisorByDivisorTester CreateGpuTester(IMersenneNumberDivisorByDivisorTester prototype) + if (tester is MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder hybridOrderPrototype) { - var gpuTester = new MersenneNumberDivisorByDivisorGpuTester - { - MinK = prototype.MinK, - }; - ((IMersenneNumberDivisorByDivisorTester)gpuTester).BatchSize = prototype.BatchSize; - return gpuTester; + RunHybridOrderFiltered(filteredPrimes, hybridOrderPrototype, workerCount, partitionSize, markComposite, clearComposite, printResult); + return; } - void ProcessPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime) + if (tester is MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder gpuOrderPrototype) + { + RunGpuOrderFiltered(filteredPrimes, gpuOrderPrototype, workerCount, partitionSize, markComposite, clearComposite, printResult); + return; + } + + if (tester is MersenneNumberDivisorByDivisorGpuTester gpuPrototype) + { + RunGpuTesterFiltered(filteredPrimes, gpuPrototype, workerCount, partitionSize, markComposite, clearComposite, printResult); + return; + } + + if (workerCount > 1) + { + throw new InvalidOperationException($"Custom by-divisor tester {tester.GetType().Name} must be thread-safe or run with --threads=1."); + } + + void ProcessCustomPrime(ulong prime) { try { Console.WriteLine($"Processing {prime}"); - IMersenneNumberDivisorByDivisorTester testerForPrime = CreateTester(tester); - string stateFile = Path.Combine(PerfectNumberConstants.ByDivisorStateDirectory, prime.ToString(CultureInfo.InvariantCulture) + ".bin"); + string stateFile = Path.Combine( + PerfectNumberConstants.ByDivisorStateDirectory, + prime.ToString(CultureInfo.InvariantCulture) + ".bin"); - BigInteger resumeK = tester.MinK; - if (File.Exists(stateFile)) - { - BigInteger lastK = File.ReadLines(stateFile) - .Where(static line => !string.IsNullOrWhiteSpace(line)) - .Select(static line => BigInteger.Parse(line, NumberStyles.None, CultureInfo.InvariantCulture)) - .Max(); + BigInteger resumeK = EnvironmentConfiguration.MinK; + if (File.Exists(stateFile)) + { + if (TryReadLastSavedK(stateFile, out BigInteger lastK)) + { + resumeK = lastK + BigInteger.One; + tester.ResumeFromState(stateFile, lastK, resumeK); + } + else + { + tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } + } + else + { + tester.ResumeFromState(stateFile, BigInteger.Zero, resumeK); + } - resumeK = lastK + BigInteger.One; - testerForPrime.ResumeFromState(lastK); - } - else - { - testerForPrime.ResumeFromState(BigInteger.Zero); - } + tester.ResetStateTracking(); - testerForPrime.MinK = resumeK; - testerForPrime.StateFilePath = stateFile; - testerForPrime.ResetStateTracking(); - - bool isPrime = testerForPrime.IsPrime(gpu, prime, out bool divisorsExhausted, out BigInteger divisor); + bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted, out BigInteger divisor); if (!isPrime) { @@ -322,70 +670,598 @@ public static class MersenneNumberDivisorByDivisorTester } } + DeterministicRandomCpu.Initialize(); + foreach (ulong prime in filteredPrimes) + { + ProcessCustomPrime(prime); + } + } + + private delegate void RunPreparedCandidates( + List filteredPrimes, + in TTester prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + where TTester : struct; + + private static void RunStructTester( + List candidates, + ref TTester tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask, + RunPreparedCandidates runPrepared) + where TTester : struct, IMersenneNumberDivisorByDivisorTester + { + List filteredPrimes = PrepareCandidatesAndGetFilteredPrimes( + candidates, + ref tester, + previousResults, + startPrime, + threadCount, + primesPerTask, + out int workerCount, + out int partitionSize); + + if (filteredPrimes.Count == 0) + { + return; + } + + TTester prototype = tester; + runPrepared(filteredPrimes, prototype, workerCount, partitionSize, markComposite, clearComposite, printResult); + } + + private static List PrepareCandidatesAndGetFilteredPrimes( + List candidates, + ref TTester tester, + Dictionary? previousResults, + ulong startPrime, + int threadCount, + int primesPerTask, + out int workerCount, + out int partitionSize) + where TTester : struct, IMersenneNumberDivisorByDivisorTester + { + if (candidates.Count == 0) + { + workerCount = 0; + partitionSize = 0; + return []; + } + + candidates.Sort(); + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + int candidateCount = candidateSpan.Length; + + bool applyStartPrime = startPrime > 0UL; + int skippedByPreviousResults = 0; + + ArrayPool pool = ThreadStaticPools.UlongPool; + if (previousResults is not null && previousResults.Count > 0) + { + int recordedCount = previousResults.Count; + ulong[] recordedCandidates = pool.Rent(recordedCount); + + Span recordedSpan = recordedCandidates.AsSpan(0, recordedCount); + int recordedIndex = 0; + foreach (ulong value in previousResults.Keys) + { + recordedSpan[recordedIndex++] = value; + } + + recordedSpan.Sort(); + + int writeIndex = 0; + int skipIndex = 0; + + for (int readIndex = 0; readIndex < candidateCount; readIndex++) + { + ulong candidate = candidateSpan[readIndex]; + + while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate) + { + skipIndex++; + } + + if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate) + { + skippedByPreviousResults++; + skipIndex++; + continue; + } + + candidateSpan[writeIndex++] = candidate; + } + + if (writeIndex < candidateCount) + { + candidates.RemoveRange(writeIndex, candidateCount - writeIndex); + candidateSpan = CollectionsMarshal.AsSpan(candidates); + candidateCount = candidateSpan.Length; + } + + pool.Return(recordedCandidates, clearArray: false); + } + + if (skippedByPreviousResults > 0) + { + Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results."); + } + + int startIndex = 0; + if (applyStartPrime) + { + startIndex = candidates.BinarySearch(startPrime); + if (startIndex < 0) + { + startIndex = ~startIndex; + } + else + { + while (startIndex > 0) + { + int previousIndex = startIndex - 1; + if (candidates[previousIndex] < startPrime) + { + break; + } + + startIndex = previousIndex; + } + } + } + + List primesToTest = candidates; + Span primesSpan = CollectionsMarshal.AsSpan(primesToTest); + ulong maxPrime = 0UL; + int primeWriteIndex = 0; + + for (int index = startIndex; index < candidateCount; index++) + { + ulong candidate = primesSpan[index]; + primesSpan[primeWriteIndex++] = candidate; + + if (candidate > maxPrime) + { + maxPrime = candidate; + } + } + + if (primeWriteIndex < primesToTest.Count) + { + primesToTest.RemoveRange(primeWriteIndex, primesToTest.Count - primeWriteIndex); + } + + if (primesToTest.Count == 0) + { + string message = applyStartPrime + ? $"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor." + : "No prime candidates remain for --mersenne=bydivisor after filtering."; + Console.WriteLine(message); + + workerCount = 0; + partitionSize = 0; + return []; + } + + int primeCount = primesToTest.Count; + ulong[] primeBatch = pool.Rent(primeCount); + ulong[] allowedMaxBatch = pool.Rent(primeCount); + List filteredPrimes = new(primeCount); + + Span primeSpan = primeBatch.AsSpan(0, primeCount); + primesToTest.CopyTo(primeSpan); + Span allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount); + + tester.PrepareCandidates(maxPrime, primeSpan, allowedMaxSpan); + + for (int i = 0; i < primeCount; i++) + { + ulong prime = primeSpan[i]; + filteredPrimes.Add(prime); + } + + pool.Return(allowedMaxBatch, clearArray: false); + pool.Return(primeBatch, clearArray: false); + + workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount; + partitionSize = primesPerTask <= 0 ? 1 : primesPerTask; + return filteredPrimes; + } + + private static void RunCpuOrderFiltered( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder prototypeCopy = prototype; + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + if (workerCount == 1) { DeterministicRandomCpu.Initialize(); - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - foreach (ulong prime in filteredPrimes) + var session = new MersenneNumberDivisorByDivisorPrimeScanSessionWithCpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + try { - ProcessPrime(gpu, prime); - } - - PrimeOrderCalculatorAccelerator.Return(gpu); - } - else - { - TaskScheduler scheduler = UnboundedTaskScheduler.Instance; - - int totalCount = filteredPrimes.Count; - int partitionSize = chunkSize < 1 ? 1 : chunkSize; - if (partitionSize > totalCount) - { - partitionSize = totalCount; - } - - int taskCount = (totalCount + partitionSize - 1) / partitionSize; - Task[] tasks = new Task[taskCount]; - var startGate = new ManualResetEventSlim(initialState: false); - int taskIndex = 0; - - for (int start = 0; start < totalCount; start += partitionSize) - { - int rangeStart = start; - int rangeEnd = rangeStart + partitionSize; - if (rangeEnd > totalCount) + foreach (ulong prime in filteredPrimes) { - rangeEnd = totalCount; + session.ProcessPrime(prime); } - - tasks[taskIndex++] = Task.Factory.StartNew( - () => - { - DeterministicRandomCpu.Initialize(); - startGate.Wait(); - GpuPrimeWorkLimiter.Acquire(); - Console.WriteLine($"Task started for range {rangeStart}"); - - PrimeOrderCalculatorAccelerator gpu = PrimeOrderCalculatorAccelerator.Rent(1); - for (int index = rangeStart; index < rangeEnd; index++) - { - ProcessPrime(gpu, filteredPrimes[index]); - } - - PrimeOrderCalculatorAccelerator.Return(gpu); - Console.WriteLine($"Task finished for range {rangeStart}"); - PrimeOrderCalculatorAccelerator.DisposeAll(); - GpuPrimeWorkLimiter.Acquire(); - }, - CancellationToken.None, - TaskCreationOptions.DenyChildAttach, - scheduler); - + } + finally + { + session.Dispose(); } - startGate.Set(); - Task.WaitAll(tasks); - startGate.Dispose(); + return; } + + TaskScheduler scheduler = UnboundedTaskScheduler.Instance; + int nextIndex = 0; + Task[] tasks = new Task[workerCount]; + var startGate = new ManualResetEventSlim(initialState: false); + + for (int taskIndex = 0; taskIndex < workerCount; taskIndex++) + { + int workerIndex = taskIndex; + tasks[taskIndex] = Task.Factory.StartNew( + () => + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorPrimeScanSessionWithCpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + startGate.Wait(); + Console.WriteLine($"Task started for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + + try + { + while (true) + { + int start = Interlocked.Add(ref nextIndex, partitionSize) - partitionSize; + if (start >= totalCount) + { + break; + } + + int end = start + partitionSize; + if (end > totalCount) + { + end = totalCount; + } + + for (int index = start; index < end; index++) + { + session.ProcessPrime(filteredPrimes[index]); + } + } + } + finally + { + session.Dispose(); + Console.WriteLine($"Task finished for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + } + }, + CancellationToken.None, + TaskCreationOptions.DenyChildAttach, + scheduler); + } + + startGate.Set(); + Task.WaitAll(tasks); + startGate.Dispose(); + } + + private static void RunHybridOrderFiltered( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder prototypeCopy = prototype; + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + + int effectiveWorkerCount = workerCount; + int gpuLimit = GpuPrimeWorkLimiter.CurrentLimit; + if (effectiveWorkerCount > gpuLimit) + { + effectiveWorkerCount = gpuLimit; + } + + if (workerCount == 1) + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder(prototypeCopy, markComposite, clearComposite, printResult); + try + { + GpuPrimeWorkLimiter.Acquire(); + foreach (ulong prime in filteredPrimes) + { + session.ProcessPrime(prime); + } + } + finally + { + GpuPrimeWorkLimiter.Release(); + session.Dispose(); + } + + return; + } + + TaskScheduler scheduler = UnboundedTaskScheduler.Instance; + int nextIndex = 0; + Task[] tasks = new Task[effectiveWorkerCount]; + var startGate = new ManualResetEventSlim(initialState: false); + + for (int taskIndex = 0; taskIndex < effectiveWorkerCount; taskIndex++) + { + int workerIndex = taskIndex; + tasks[taskIndex] = Task.Factory.StartNew( + () => + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorPrimeScanSessionWithHybridOrder(prototypeCopy, markComposite, clearComposite, printResult); + startGate.Wait(); + Console.WriteLine($"Task started for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + + try + { + GpuPrimeWorkLimiter.Acquire(); + while (true) + { + int start = Interlocked.Add(ref nextIndex, partitionSize) - partitionSize; + if (start >= totalCount) + { + break; + } + + int end = start + partitionSize; + if (end > totalCount) + { + end = totalCount; + } + + for (int index = start; index < end; index++) + { + session.ProcessPrime(filteredPrimes[index]); + } + } + } + finally + { + GpuPrimeWorkLimiter.Release(); + session.Dispose(); + Console.WriteLine($"Task finished for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + } + }, + CancellationToken.None, + TaskCreationOptions.DenyChildAttach, + scheduler); + } + + startGate.Set(); + Task.WaitAll(tasks); + startGate.Dispose(); + } + + private static void RunGpuOrderFiltered( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder prototypeCopy = prototype; + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + + int effectiveWorkerCount = workerCount; + int gpuLimit = GpuPrimeWorkLimiter.CurrentLimit; + if (effectiveWorkerCount > gpuLimit) + { + effectiveWorkerCount = gpuLimit; + } + + if (workerCount == 1) + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + try + { + GpuPrimeWorkLimiter.Acquire(); + foreach (ulong prime in filteredPrimes) + { + session.ProcessPrime(prime); + } + } + finally + { + GpuPrimeWorkLimiter.Release(); + session.Dispose(); + } + + return; + } + + TaskScheduler scheduler = UnboundedTaskScheduler.Instance; + int nextIndex = 0; + Task[] tasks = new Task[effectiveWorkerCount]; + var startGate = new ManualResetEventSlim(initialState: false); + + for (int taskIndex = 0; taskIndex < effectiveWorkerCount; taskIndex++) + { + int workerIndex = taskIndex; + tasks[taskIndex] = Task.Factory.StartNew( + () => + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorPrimeScanSessionWithGpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + startGate.Wait(); + Console.WriteLine($"Task started for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + + try + { + GpuPrimeWorkLimiter.Acquire(); + while (true) + { + int start = Interlocked.Add(ref nextIndex, partitionSize) - partitionSize; + if (start >= totalCount) + { + break; + } + + int end = start + partitionSize; + if (end > totalCount) + { + end = totalCount; + } + + for (int index = start; index < end; index++) + { + session.ProcessPrime(filteredPrimes[index]); + } + } + } + finally + { + GpuPrimeWorkLimiter.Release(); + session.Dispose(); + Console.WriteLine($"Task finished for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + } + }, + CancellationToken.None, + TaskCreationOptions.DenyChildAttach, + scheduler); + } + + startGate.Set(); + Task.WaitAll(tasks); + startGate.Dispose(); + } + + private static void RunGpuTesterFiltered( + List filteredPrimes, + MersenneNumberDivisorByDivisorGpuTester prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + + int effectiveWorkerCount = workerCount; + int gpuLimit = GpuPrimeWorkLimiter.CurrentLimit; + if (effectiveWorkerCount > gpuLimit) + { + effectiveWorkerCount = gpuLimit; + } + + if (workerCount == 1) + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorGpuPrimeScanSession(prototype, markComposite, clearComposite, printResult); + try + { + GpuPrimeWorkLimiter.Acquire(); + foreach (ulong prime in filteredPrimes) + { + session.ProcessPrime(prime); + } + } + finally + { + GpuPrimeWorkLimiter.Release(); + session.Dispose(); + } + + return; + } + + TaskScheduler scheduler = UnboundedTaskScheduler.Instance; + int nextIndex = 0; + Task[] tasks = new Task[effectiveWorkerCount]; + var startGate = new ManualResetEventSlim(initialState: false); + + for (int taskIndex = 0; taskIndex < effectiveWorkerCount; taskIndex++) + { + int workerIndex = taskIndex; + tasks[taskIndex] = Task.Factory.StartNew( + () => + { + DeterministicRandomCpu.Initialize(); + var session = new MersenneNumberDivisorByDivisorGpuPrimeScanSession(prototype, markComposite, clearComposite, printResult); + startGate.Wait(); + Console.WriteLine($"Task started for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + + try + { + GpuPrimeWorkLimiter.Acquire(); + while (true) + { + int start = Interlocked.Add(ref nextIndex, partitionSize) - partitionSize; + if (start >= totalCount) + { + break; + } + + int end = start + partitionSize; + if (end > totalCount) + { + end = totalCount; + } + + for (int index = start; index < end; index++) + { + session.ProcessPrime(filteredPrimes[index]); + } + } + } + finally + { + GpuPrimeWorkLimiter.Release(); + session.Dispose(); + Console.WriteLine($"Task finished for worker {workerIndex.ToString(CultureInfo.InvariantCulture)}"); + } + }, + CancellationToken.None, + TaskCreationOptions.DenyChildAttach, + scheduler); + } + + startGate.Set(); + Task.WaitAll(tasks); + startGate.Dispose(); } } diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 6e4d693..45f619d 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -12,14 +12,14 @@ public static class PerfectNumberConstants public const int DefaultThreadPoolCapacity = 30_000; public const int DefaultSpecialMaxFactorCapacity = 1024; public const ulong ExtraDivisorCycleSearchLimit = 64UL; - public static int GpuRatio = 16384 - 1; + public const int GpuRatio = 16384 - 1; public const int MaxQForDivisorCycles = 4_000_000; - public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; + public static readonly uint PrimesLimit = 1_000_000; public const int DefaultStringBuilderCapacity = 10_240; public const int MaxOddPowersCount = 128; public const int PooledArrayThreshold = 64; public const int Pow2WindowSize = 8; - public static int RollingAccelerators = 298; //SharedGpuContext.Device.MaxNumThreadsPerGroup; + public const int RollingAccelerators = 298; public const int ThreadsByAccelerator = 1; public const int ByDivisorStateSaveInterval = 200_000; public const string ByDivisorStateDirectory = "Checks"; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs index 3339b62..e0b936e 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs @@ -328,12 +328,13 @@ internal static partial class PrimeOrderCalculator [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static bool RunOnCpu() { - if (PerfectNumberConstants.GpuRatio == 0) + int gpuRatio = EnvironmentConfiguration.GpuRatio; + if (gpuRatio == 0) { return false; } - if (PerfectNumberConstants.GpuRatio == 1) + if (gpuRatio == 1) { while (true) { @@ -352,16 +353,16 @@ internal static partial class PrimeOrderCalculator } int cpuCount = Interlocked.Increment(ref _cpuCount); - if (cpuCount == PerfectNumberConstants.GpuRatio) + if (cpuCount == gpuRatio) { - Interlocked.Add(ref _cpuCount, -PerfectNumberConstants.GpuRatio); + Interlocked.Add(ref _cpuCount, -gpuRatio); } - else if (cpuCount > PerfectNumberConstants.GpuRatio) + else if (cpuCount > gpuRatio) { - cpuCount -= PerfectNumberConstants.GpuRatio; + cpuCount -= gpuRatio; } - return cpuCount != PerfectNumberConstants.GpuRatio; + return cpuCount != gpuRatio; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] diff --git a/PerfectNumbers.Core/TaskThreadPool.cs b/PerfectNumbers.Core/TaskThreadPool.cs index f746ee4..2f131b3 100644 --- a/PerfectNumbers.Core/TaskThreadPool.cs +++ b/PerfectNumbers.Core/TaskThreadPool.cs @@ -153,9 +153,9 @@ internal sealed class TaskThreadPool return false; } - workItem = _pendingTaskQueue.Pop(); - if (workItem is { }) + if (_pendingTaskQueue.TryPop(out WorkItem dequeued)) { + workItem = dequeued; Interlocked.Decrement(ref _pendingTaskCount); return true; } diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 61a2dad..c6af7d6 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -127,24 +127,24 @@ namespace PerfectNumbers.Core get => _gpuUInt128Pool ??= ArrayPool.Create(); } - [ThreadStatic] - private static MersenneCpuDivisorScanSession? _mersenneCpuDivisorSession; + // [ThreadStatic] + // private static MersenneCpuDivisorScanSession? _mersenneCpuDivisorSession; - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - internal static MersenneCpuDivisorScanSession RentMersenneCpuDivisorSession(Gpu.Accelerators.PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) - { - if (_mersenneCpuDivisorSession is { } session) - { - session.Configure(gpu, orderDevice); - _mersenneCpuDivisorSession = null; - return session; - } + // [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + // internal static MersenneCpuDivisorScanSession RentMersenneCpuDivisorSession(Gpu.Accelerators.PrimeOrderCalculatorAccelerator gpu, ComputationDevice orderDevice) + // { + // if (_mersenneCpuDivisorSession is { } session) + // { + // session.Configure(gpu, orderDevice); + // _mersenneCpuDivisorSession = null; + // return session; + // } - return new MersenneCpuDivisorScanSession(gpu, orderDevice); - } + // return new MersenneCpuDivisorScanSession(gpu, orderDevice); + // } - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - internal static void ReturnMersenneCpuDivisorSession(MersenneCpuDivisorScanSession session) => _mersenneCpuDivisorSession = session; + // [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + // internal static void ReturnMersenneCpuDivisorSession(MersenneCpuDivisorScanSession session) => _mersenneCpuDivisorSession = session; [ThreadStatic] private static ArrayPool? _uint128Pool;