From ed299edf56cc646306d01a2d5e9140f17052a38a Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 18 Dec 2025 20:10:20 +0100 Subject: [PATCH] MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder, MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder & MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder converted to device dependent template MersenneNumberDivisorByDivisorCpuTesterWithTemplate. --- ...erDivisorByDivisorCpuTesterWithCpuOrder.cs | 585 ----------------- ...erDivisorByDivisorCpuTesterWithGpuOrder.cs | 600 ------------------ ...rDivisorByDivisorCpuTesterWithTemplate.cs} | 166 ++++- 3 files changed, 143 insertions(+), 1208 deletions(-) delete mode 100644 PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs delete mode 100644 PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder.cs rename PerfectNumbers.Core/Cpu/{MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder.cs => MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs} (78%) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs deleted file mode 100644 index 54c0a26..0000000 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder.cs +++ /dev/null @@ -1,585 +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 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 deleted file mode 100644 index 937e3ac..0000000 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder.cs +++ /dev/null @@ -1,600 +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 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/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs similarity index 78% rename from PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder.cs rename to PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs index 28f0beb..9b33887 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs @@ -6,11 +6,15 @@ using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core.Cpu; -public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersenneNumberDivisorByDivisorTester +[DeviceDependentTemplate(typeof(ComputationDevice), suffix: "Order")] +public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneNumberDivisorByDivisorTester { public BigInteger DivisorLimit; public BigInteger MinK = EnvironmentConfiguration.MinK; + +#if DEVICE_GPU || DEVICE_HYBRID public readonly PrimeOrderCalculatorAccelerator Accelerator = PrimeOrderCalculatorAccelerator.Rent(1); +#endif #pragma warning disable CS8618 // StateFilePath is always set on EvenPerfectBitScanner execution path when ResumeFromState is called. private string StateFilePath; @@ -21,11 +25,13 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen 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; @@ -81,9 +87,17 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public readonly IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) { +#if DEVICE_GPU + var session = new MersenneCpuDivisorScanSessionWithGpuOrder(); + session.Configure(gpu); + return session; +#elif DEVICE_HYBRID var session = new MersenneCpuDivisorScanSessionWithHybridOrder(); session.Configure(gpu); return session; +#else + return new MersenneCpuDivisorScanSessionWithCpuOrder(); +#endif } private bool CheckDivisors( @@ -254,32 +268,69 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen return composite; } +#if DEVICE_GPU || DEVICE_HYBRID PrimeOrderCalculatorAccelerator gpu = Accelerator; +#endif var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); + ulong computedCycle; + bool primeOrderFailed; + ulong candidate; + ulong divisorCycle; + bool computed; while (divisor.CompareTo(limit) <= 0) { if (residueStepper.IsAdmissible()) { - ulong candidate = divisor.Low; + 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) +#if DEVICE_GPU + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( + candidate, + gpu, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#elif DEVICE_HYBRID + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( + candidate, + gpu, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#else + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( + candidate, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#endif + if (!computed || 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. +#if DEVICE_GPU + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( + candidate, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); +#elif DEVICE_HYBRID divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( candidate, gpu, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); +#else + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( + candidate, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); +#endif } else { @@ -312,6 +363,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen return false; } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisorsLarge( ulong prime, BigInteger allowedMax, @@ -402,6 +454,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen return false; } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisors64Range( ulong prime, ulong step, @@ -452,6 +505,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen out foundDivisor); } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisors64( ulong prime, ulong step, @@ -462,27 +516,63 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen out bool processedAll, out ulong foundDivisor) { +#if DEVICE_GPU || DEVICE_HYBRID PrimeOrderCalculatorAccelerator gpu = Accelerator; +#endif + ulong computedCycle; + bool primeOrderFailed; bool canAdvance = step <= limit; + ulong divisorCycle; + bool computed; if (!canAdvance) { if (divisor <= limit && residueStepper.IsAdmissible()) { MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - ulong divisorCycle; - if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( +#if DEVICE_GPU + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( + divisor, + gpu, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#elif DEVICE_HYBRID + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( + divisor, + gpu, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#else + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( + divisor, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#endif + if (!computed || computedCycle == 0UL) + { +#if DEVICE_GPU + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( divisor, gpu, - prime, divisorData, - out ulong computedCycle, - out bool primeOrderFailed) || computedCycle == 0UL) - { + skipPrimeOrderHeuristic: primeOrderFailed); +#elif DEVICE_HYBRID divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( divisor, gpu, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); +#else + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( + divisor, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); +#endif } else { @@ -520,20 +610,50 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen if (residueStepper.IsAdmissible()) { MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - ulong divisorCycle; - if (!MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( +#if DEVICE_GPU + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( + divisor, + gpu, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#elif DEVICE_HYBRID + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( + divisor, + gpu, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#else + computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( + divisor, + prime, + divisorData, + out computedCycle, + out primeOrderFailed); +#endif + if (!computed || computedCycle == 0UL) + { +#if DEVICE_GPU + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( divisor, gpu, - prime, divisorData, - out ulong computedCycle, - out bool primeOrderFailed) || computedCycle == 0UL) - { + skipPrimeOrderHeuristic: primeOrderFailed); +#elif DEVICE_HYBRID divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( divisor, gpu, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); +#else + divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( + divisor, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); +#endif } else { @@ -568,10 +688,10 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder() : IMersen } } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private void RecordState(BigInteger k) { string path = StateFilePath; - if (k <= _lastSavedK) { return;