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 : IMersenneNumberDivisorByDivisorTester { private BigInteger _divisorLimit; private int _batchSize = 1_024; private BigInteger _minK = BigInteger.One; private readonly ComputationDevice _orderDevice; private readonly TryCycleLengthDelegate _tryCalculateCycleLengthForExponent; private readonly CalculateCycleLengthDelegate _calculateCycleLength; private string? _stateFilePath; private int _stateCounter; private BigInteger _lastSavedK; public ComputationDevice OrderDevice => _orderDevice; 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 MersenneNumberDivisorByDivisorCpuTester(ComputationDevice orderDevice = ComputationDevice.Hybrid) { _orderDevice = orderDevice; _tryCalculateCycleLengthForExponent = orderDevice switch { ComputationDevice.Gpu => TryCalculateCycleLengthForExponentGpu, ComputationDevice.Hybrid => TryCalculateCycleLengthForExponentHybrid, _ => TryCalculateCycleLengthForExponentCpu, }; _calculateCycleLength = orderDevice switch { ComputationDevice.Gpu => CalculateCycleLengthGpu, ComputationDevice.Hybrid => CalculateCycleLengthHybrid, _ => CalculateCycleLengthCpu, }; } public int BatchSize { 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 processedAll; bool composite = CheckDivisors( gpu, prime, allowedMax, _minK, out 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); var divisorPool = MontgomeryDivisorDataPool.Shared; while (divisor.CompareTo(limit) <= 0) { if (residueStepper.IsAdmissible()) { ulong candidate = divisor.Low; MontgomeryDivisorData divisorData = divisorPool.FromModulus(candidate); ulong divisorCycle; // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. if (!_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; } divisorPool.Return(divisorData); 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; FixedCapacityStack divisorPool = MontgomeryDivisorDataPool.Shared; if (!canAdvance) { if (divisor <= limit && residueStepper.IsAdmissible()) { MontgomeryDivisorData divisorData = divisorPool.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; } divisorPool.Return(divisorData); 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 = divisorPool.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; } divisorPool.Return(divisorData); 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; } }