using System; using System.Numerics; using System.Runtime.CompilerServices; using System.Globalization; using System.Collections.Generic; using System.Collections.Concurrent; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core.Cpu; public enum DivisorSet { OneByOne, Pow2Groups, } [EnumDependentTemplate(typeof(DivisorSet))] [DeviceDependentTemplate(typeof(ComputationDevice))] [NameSuffix("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; #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 bool _pow2MinusOneChecked; private string _pow2MinusOneStateFilePath = string.Empty; private GpuUInt128WorkSet _divisorScanGpuWorkSet; private static readonly ConcurrentDictionary _checkedPow2MinusOne = new(); #if DivisorSet_Pow2Groups private const int PercentScale = DivisorByDivisorConfig.PercentScale; private static readonly int[] SpecialPercentPromille = DivisorByDivisorConfig.SpecialPercentPromille; private static readonly (int Start, int End)[] PercentGroupsPromille = DivisorByDivisorConfig.PercentGroupsPromille; private static long _specialHits; private static long _groupHits; private enum Pow2Phase { None, Special, Groups, } private Pow2Phase _currentPow2Phase; private string _specialStateFilePath = string.Empty; private string _groupsStateFilePath = string.Empty; private int _specialStateCounter; private int _groupsStateCounter; private BigInteger _specialLastSavedK; private BigInteger _groupsLastSavedK; private BigInteger _specialResumeK; private BigInteger _groupsResumeK; private bool _preparedSpecialsInitialized; private int _preparedSpecialRangePromille; private int[] _preparedSpecialSingles = Array.Empty(); private (int Start, int End)[] _preparedSpecialRanges = Array.Empty<(int Start, int End)>(); #endif [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public void ResetStateTracking() { _stateCounter = 0; _pow2MinusOneChecked = false; #if DivisorSet_Pow2Groups _specialStateCounter = 0; _groupsStateCounter = 0; _currentPow2Phase = Pow2Phase.None; _preparedSpecialsInitialized = false; #endif } [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; _pow2MinusOneStateFilePath = stateFile + ".pow2minus1"; _pow2MinusOneChecked = File.Exists(_pow2MinusOneStateFilePath); #if DivisorSet_Pow2Groups _specialStateFilePath = stateFile + ".special"; _groupsStateFilePath = stateFile + ".groups"; _specialStateCounter = 0; _groupsStateCounter = 0; if (!MersenneNumberDivisorByDivisorTester.TryReadLastSavedK(_specialStateFilePath, out _specialLastSavedK)) { _specialLastSavedK = lastSavedK; } if (!MersenneNumberDivisorByDivisorTester.TryReadLastSavedK(_groupsStateFilePath, out _groupsLastSavedK)) { _groupsLastSavedK = lastSavedK; } _specialResumeK = _specialLastSavedK + BigInteger.One; _groupsResumeK = _groupsLastSavedK + BigInteger.One; _currentPow2Phase = Pow2Phase.None; _preparedSpecialsInitialized = false; #endif } [MethodImpl(MethodImplOptions.AggressiveOptimization)] public bool IsPrime(ulong prime, out bool divisorsExhausted, out BigInteger divisor) { BigInteger allowedMax = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisorBig(prime, DivisorLimit); if (TryFindPowerOfTwoMinusOneDivisor(prime, allowedMax, out divisor)) { divisorsExhausted = true; return false; } // 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 divisorsExhausted, out divisor); if (composite) { divisorsExhausted = true; return false; } 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) { #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 } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisors( ulong prime, BigInteger allowedMax, BigInteger minK, out bool processedAll, out BigInteger foundDivisor) { // TODO: Is minK ever < 1? Can this be optimized - normalizedMinK is always >= 1 after this BigInteger normalizedMinK = minK >= BigInteger.One ? minK : BigInteger.One, step = ((BigInteger)prime) << 1, firstDivisor = (step * normalizedMinK) + BigInteger.One; bool primeLastIsSeven = (prime & 3UL) == 3UL; ushort primeDecimalMask = DivisorGenerator.GetDecimalMask(primeLastIsSeven); bool fits64 = normalizedMinK <= ulong.MaxValue && firstDivisor <= ulong.MaxValue, unlimited = allowedMax.IsZero, // This checks if the limit fits 64-bit composite64 = unlimited || allowedMax <= ulong.MaxValue; ulong allowedMax64, foundDivisor64; #if DivisorSet_Pow2Groups bool compositePow2 = CheckDivisorsPow2Groups( prime, primeDecimalMask, allowedMax, normalizedMinK, step, out processedAll, out foundDivisor); return compositePow2; #else if (fits64 && composite64) { allowedMax64 = unlimited ? ulong.MaxValue : (ulong)allowedMax; composite64 = CheckDivisors64Bit( prime, primeDecimalMask, allowedMax64, (ulong)normalizedMinK, out processedAll, out foundDivisor64); foundDivisor = foundDivisor64; return composite64; } if (fits64 && allowedMax > ulong.MaxValue) { allowedMax64 = ulong.MaxValue; composite64 = CheckDivisors64Bit( prime, primeDecimalMask, allowedMax64, (ulong)normalizedMinK, out processedAll, out foundDivisor64); if (composite64) { foundDivisor = foundDivisor64; return true; } BigInteger iterations = ((BigInteger)(allowedMax64 - (ulong)firstDivisor) / (ulong)step) + BigInteger.One; BigInteger nextK = normalizedMinK + iterations; return CheckDivisorsLarge( prime, primeDecimalMask, allowedMax, nextK, step, out processedAll, out foundDivisor); } #endif return CheckDivisorsLarge( prime, primeDecimalMask, allowedMax, normalizedMinK, step, out processedAll, out foundDivisor); } #if DivisorSet_Pow2Groups private static readonly (int, int)[] _emptyPow2Groups = Array.Empty<(int, int)>(); [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisorsPow2Groups( ulong prime, ushort primeDecimalMask, BigInteger allowedMax, BigInteger minK, BigInteger step, out bool processedAll, out BigInteger foundDivisor) { foundDivisor = BigInteger.Zero; processedAll = true; if (allowedMax.IsZero) { processedAll = false; return false; } BigInteger normalizedMinK = minK >= BigInteger.One ? minK : BigInteger.One; BigInteger maxKAllowed = (allowedMax - BigInteger.One) / step; if (maxKAllowed < normalizedMinK) { return false; } ReadOnlySpan specialPercents = SpecialPercentPromille; ReadOnlySpan<(int Start, int End)> groups = PercentGroupsPromille; // Segment 1: ulong range if (normalizedMinK <= BigIntegerNumbers.UlongMaxValue) { ulong startK64 = (ulong)normalizedMinK; ulong maxK64 = maxKAllowed > BigIntegerNumbers.UlongMaxValue ? ulong.MaxValue : (ulong)maxKAllowed; ulong allowedMax64 = allowedMax > BigIntegerNumbers.UlongMaxValue ? ulong.MaxValue : (ulong)allowedMax; bool specialFound = ProcessPow2Phase( prime, allowedMax64, prime << 1, _specialResumeK > BigInteger.Zero ? _specialResumeK : startK64, startK64, Pow2Phase.Special, _emptyPow2Groups, specialPercents, ref processedAll, out ulong foundDiv64); if (specialFound) { foundDivisor = foundDiv64; return true; } bool groupFound = ProcessPow2Phase( prime, allowedMax64, prime << 1, _groupsResumeK > BigInteger.Zero ? _groupsResumeK : startK64, startK64, Pow2Phase.Groups, groups, ReadOnlySpan.Empty, ref processedAll, out foundDiv64); if (groupFound) { foundDivisor = foundDiv64; return true; } } // Segment 2: UInt128 range BigInteger start128 = BigInteger.Max(normalizedMinK, BigIntegerNumbers.UlongMaxValuePlusOne); if (start128 <= maxKAllowed && start128 <= BigIntegerNumbers.UInt128MaxValue) { BigInteger end128 = BigInteger.Min(maxKAllowed, BigIntegerNumbers.UInt128MaxValue); bool specialFound128 = ProcessPow2PhaseBig( prime, primeDecimalMask, allowedMax, step, _specialResumeK > start128 ? _specialResumeK : start128, normalizedMinK, start128, end128, Pow2Phase.Special, _emptyPow2Groups, specialPercents, ref processedAll, out BigInteger foundBig); if (specialFound128) { foundDivisor = foundBig; return true; } bool groupFound128 = ProcessPow2PhaseBig( prime, primeDecimalMask, allowedMax, step, _groupsResumeK > start128 ? _groupsResumeK : start128, normalizedMinK, start128, end128, Pow2Phase.Groups, groups, ReadOnlySpan.Empty, ref processedAll, out foundBig); if (groupFound128) { foundDivisor = foundBig; return true; } } // Segment 3: BigInteger range BigInteger startBig = BigInteger.Max(normalizedMinK, BigIntegerNumbers.UInt128MaxValuePlusOne); if (startBig <= maxKAllowed) { bool specialFoundBig = ProcessPow2PhaseBig( prime, primeDecimalMask, allowedMax, step, _specialResumeK > startBig ? _specialResumeK : startBig, normalizedMinK, startBig, maxKAllowed, Pow2Phase.Special, _emptyPow2Groups, specialPercents, ref processedAll, out BigInteger foundBig); if (specialFoundBig) { foundDivisor = foundBig; return true; } bool groupFoundBig = ProcessPow2PhaseBig( prime, primeDecimalMask, allowedMax, step, _groupsResumeK > startBig ? _groupsResumeK : startBig, normalizedMinK, startBig, maxKAllowed, Pow2Phase.Groups, groups, ReadOnlySpan.Empty, ref processedAll, out foundBig); if (groupFoundBig) { foundDivisor = foundBig; return true; } } return false; } #endif [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisors64Bit( ulong prime, ushort primeDecimalMask, ulong allowedMax, ulong minK, out bool processedAll, out ulong foundDivisor) { foundDivisor = 0UL; #if DivisorSet_Pow2Groups ulong step64 = prime << 1; if (step64 == 0UL) { processedAll = true; return false; } processedAll = true; BigInteger stepBig = step64; BigInteger maxKAllowed = allowedMax == 0UL ? BigInteger.Zero : (((BigInteger)allowedMax - BigInteger.One) / stepBig); ReadOnlySpan specialPercents = SpecialPercentPromille; ReadOnlySpan<(int Start, int End)> groups = PercentGroupsPromille; // Segment 1: ulong range ulong max64K = maxKAllowed <= ulong.MaxValue ? (ulong)maxKAllowed : ulong.MaxValue; if (minK <= max64K) { ulong max64 = allowedMax == 0UL ? ulong.MaxValue : allowedMax; bool specialFound = ProcessPow2Phase( prime, max64, step64, _specialResumeK > BigInteger.Zero ? _specialResumeK : minK, minK, Pow2Phase.Special, _emptyPow2Groups, specialPercents, ref processedAll, out foundDivisor); if (specialFound) { return true; } bool groupFound = ProcessPow2Phase( prime, max64, step64, _groupsResumeK > BigInteger.Zero ? _groupsResumeK : minK, minK, Pow2Phase.Groups, groups, ReadOnlySpan.Empty, ref processedAll, out foundDivisor); if (groupFound) { return true; } } // Segment 2: UInt128 range BigInteger start128 = BigInteger.Max(minK, BigIntegerNumbers.UlongMaxValuePlusOne); if (maxKAllowed > BigIntegerNumbers.UlongMaxValue && start128 <= maxKAllowed) { BigInteger end128 = BigInteger.Min(maxKAllowed, BigIntegerNumbers.UInt128MaxValue); bool specialFound128 = ProcessPow2PhaseBig( prime, primeDecimalMask, (BigInteger)allowedMax == BigInteger.Zero ? BigInteger.Zero : (BigInteger)allowedMax, stepBig, _specialResumeK > start128 ? _specialResumeK : start128, minK, start128, end128, Pow2Phase.Special, _emptyPow2Groups, specialPercents, ref processedAll, out BigInteger foundBig); if (specialFound128) { foundDivisor = foundBig > ulong.MaxValue ? ulong.MaxValue : (ulong)foundBig; return true; } bool groupFound128 = ProcessPow2PhaseBig( prime, primeDecimalMask, (BigInteger)allowedMax == BigInteger.Zero ? BigInteger.Zero : (BigInteger)allowedMax, stepBig, _groupsResumeK > start128 ? _groupsResumeK : start128, minK, start128, end128, Pow2Phase.Groups, groups, ReadOnlySpan.Empty, ref processedAll, out foundBig); if (groupFound128) { foundDivisor = foundBig > ulong.MaxValue ? ulong.MaxValue : (ulong)foundBig; return true; } } // Segment 3: BigInteger range beyond UInt128 BigInteger startBig = BigInteger.Max(minK, BigIntegerNumbers.UInt128MaxValuePlusOne); if (maxKAllowed > BigIntegerNumbers.UInt128MaxValue && startBig <= maxKAllowed) { bool specialFoundBig = ProcessPow2PhaseBig( prime, primeDecimalMask, (BigInteger)allowedMax == BigInteger.Zero ? BigInteger.Zero : (BigInteger)allowedMax, stepBig, _specialResumeK > startBig ? _specialResumeK : startBig, minK, startBig, maxKAllowed, Pow2Phase.Special, _emptyPow2Groups, specialPercents, ref processedAll, out BigInteger foundBig); if (specialFoundBig) { foundDivisor = foundBig > ulong.MaxValue ? ulong.MaxValue : (ulong)foundBig; return true; } bool groupFoundBig = ProcessPow2PhaseBig( prime, primeDecimalMask, (BigInteger)allowedMax == BigInteger.Zero ? BigInteger.Zero : (BigInteger)allowedMax, stepBig, _groupsResumeK > startBig ? _groupsResumeK : startBig, minK, startBig, maxKAllowed, Pow2Phase.Groups, groups, ReadOnlySpan.Empty, ref processedAll, out foundBig); if (groupFoundBig) { foundDivisor = foundBig > ulong.MaxValue ? ulong.MaxValue : (ulong)foundBig; return true; } } return false; #else // 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 currentK = minK, // We're reusing candidate variable as stepLow to limit registry pressure candidate = step.Low, // We're reusing computedCycle variable as lowerEnd to limit registry pressure computedCycle, // We're reusing divisorCycle variable as maxK to limit registry pressure divisorCycle; // We're reusing computed variable as processedTop to limit registry pressure bool computed, // We're reusing primeOrderFailed variable as processedBottom to limit registry pressure primeOrderFailed, composite; if (step.High == 0UL && limit.High == 0UL) { divisorCycle = allowedMax > 0UL ? (allowedMax - 1UL) / candidate : 0UL; if (divisorCycle == 0UL) { processedAll = true; foundDivisor = 0UL; return false; } computed = true; primeOrderFailed = true; composite = false; if (minK <= divisorCycle) { composite = CheckDivisors64Range( prime, candidate, allowedMax, minK, divisorCycle, ref currentK, out computed, out foundDivisor); if (composite) { processedAll = true; return true; } } computedCycle = minK > 1UL ? minK - 1UL : 0UL; if (computedCycle >= 1UL && minK == 1UL) { computedCycle = Math.Min(computedCycle, divisorCycle); currentK = 1UL; composite = CheckDivisors64Range( prime, candidate, allowedMax, 1UL, computedCycle, ref currentK, out primeOrderFailed, out foundDivisor); } processedAll = computed && primeOrderFailed; return composite; } #if DEVICE_GPU || DEVICE_HYBRID PrimeOrderCalculatorAccelerator gpu = Accelerator; #endif var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); while (divisor.CompareTo(limit) <= 0) { if (residueStepper.IsAdmissible()) { candidate = divisor.Low; MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); // 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 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) { divisorCycle = computedCycle; } else { // 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 } 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; #endif } #if DivisorSet_Pow2Groups [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool ProcessPow2Phase( ulong prime, ulong max64, ulong step64, BigInteger resumeK, ulong minK, Pow2Phase phase, ReadOnlySpan<(int Start, int End)> ranges, ReadOnlySpan singles, ref bool processedAll, out ulong foundDivisor) { foundDivisor = 0UL; int specialRangePromille = EnvironmentConfiguration.ByDivisorSpecialRange; if (specialRangePromille < 0) { specialRangePromille = 0; } specialRangePromille *= 100; // convert percent to promille-of-percentScale (10_000 = 100%) EnsurePreparedSpecials(specialRangePromille); if (resumeK > ulong.MaxValue) { processedAll = false; return false; } BigInteger effectiveResume = resumeK < minK ? minK : resumeK; BigInteger firstDivisorBig = (((BigInteger)prime) << 1) * effectiveResume + BigInteger.One; if (firstDivisorBig > max64) { processedAll = false; return false; } ulong firstDivisor = (ulong)firstDivisorBig; int startBit = BitOperations.Log2(firstDivisor); int maxBit = BitOperations.Log2(max64); ulong currentK = (ulong)effectiveResume; _currentPow2Phase = phase; FlushPow2State(effectiveResume, phase); ReadOnlySpan effectiveSingles = singles; ReadOnlySpan<(int Start, int End)> effectiveRanges = ranges; if (phase == Pow2Phase.Special) { effectiveSingles = _preparedSpecialSingles; effectiveRanges = _preparedSpecialRanges; } try { for (int bit = startBit; bit <= maxBit; bit++) { ulong lower = 1UL << bit; ulong upperExclusive = bit == 63 ? ulong.MaxValue : lower << 1; if (lower > max64) { break; } BigInteger lowerBig = lower; for (int i = 0; i < effectiveSingles.Length; i++) { int startProm = effectiveSingles[i]; int endProm = effectiveSingles[i]; if (specialRangePromille > 0) { startProm = Math.Max(0, effectiveSingles[i] - specialRangePromille); endProm = Math.Min(PercentScale, effectiveSingles[i] + specialRangePromille); } if (specialRangePromille == 0) { int adjustedEndProm = Math.Min(PercentScale, startProm + 1); BigInteger scaledStart = (PercentScale + startProm) * lowerBig; BigInteger scaledEnd = (PercentScale + adjustedEndProm) * lowerBig; BigInteger qMin = (scaledStart + (PercentScale - 1)) / PercentScale; BigInteger qMax = scaledEnd / PercentScale; if ((qMin & 1) == 0) { qMin += 1; } if ((qMax & 1) == 0) { qMax -= 1; } if (qMin < 3 || qMin > qMax) { continue; } if (qMin < lowerBig) { qMin = lowerBig | 1; } if (qMax >= upperExclusive) { qMax = upperExclusive - 1; if ((qMax & 1) == 0) { qMax -= 1; } } if (qMin > qMax || qMin > max64) { continue; } if (qMax > max64) { qMax = max64; if ((qMax & 1) == 0) { qMax -= 1; } } BigInteger stepBig = step64; BigInteger kMinBig = (qMin - 1 + (stepBig - 1)) / stepBig; BigInteger kMaxBig = (qMax - 1) / stepBig; if (kMinBig < minK || kMinBig < effectiveResume) { kMinBig = BigInteger.Max(BigInteger.Max(minK, effectiveResume), kMinBig); } if (kMinBig > kMaxBig) { continue; } if (kMaxBig > ulong.MaxValue) { kMaxBig = ulong.MaxValue; } ulong kStart = (ulong)kMinBig; ulong kEnd = (ulong)kMaxBig; bool composite = CheckDivisors64Range( prime, step64, max64, kStart, kEnd, ref currentK, out bool processedRange, out foundDivisor); if (composite) { RecordPow2Hit(phase); FlushPow2State(currentK, phase); _currentPow2Phase = Pow2Phase.None; return true; } if (!processedRange) { processedAll = false; } } else { BigInteger scaledStart = (PercentScale + startProm) * lowerBig; BigInteger scaledEnd = (PercentScale + endProm) * lowerBig; BigInteger qMin = (scaledStart + (PercentScale - 1)) / PercentScale; BigInteger qMax = scaledEnd / PercentScale; if ((qMin & 1) == 0) { qMin += 1; } if ((qMax & 1) == 0) { qMax -= 1; } if (qMin < 3 || qMin > qMax) { continue; } if (qMin < lowerBig) { qMin = lowerBig | 1; } if (qMax >= upperExclusive) { qMax = upperExclusive - 1; if ((qMax & 1) == 0) { qMax -= 1; } } if (qMin > qMax || qMin > max64) { continue; } if (qMax > max64) { qMax = max64; if ((qMax & 1) == 0) { qMax -= 1; } } BigInteger stepBig = step64; BigInteger kMinBig = (qMin - 1 + (stepBig - 1)) / stepBig; BigInteger kMaxBig = (qMax - 1) / stepBig; if (kMinBig < minK || kMinBig < effectiveResume) { kMinBig = BigInteger.Max(minK, effectiveResume); } if (kMinBig > kMaxBig) { continue; } if (kMaxBig > ulong.MaxValue) { kMaxBig = ulong.MaxValue; } ulong kStart = (ulong)kMinBig; ulong kEnd = (ulong)kMaxBig; bool composite = CheckDivisors64Range( prime, step64, max64, kStart, kEnd, ref currentK, out bool processedRange, out foundDivisor); if (composite) { RecordPow2Hit(phase); FlushPow2State(currentK, phase); _currentPow2Phase = Pow2Phase.None; return true; } if (!processedRange) { processedAll = false; } } } for (int i = 0; i < effectiveRanges.Length; i++) { int start = effectiveRanges[i].Start; int end = effectiveRanges[i].End; BigInteger scaledStart = (PercentScale + start) * lowerBig; BigInteger scaledEnd = (PercentScale + end) * lowerBig; BigInteger qMin = (scaledStart + (PercentScale - 1)) / PercentScale; BigInteger qMax = scaledEnd / PercentScale; if ((qMin & 1) == 0) { qMin += 1; } if ((qMax & 1) == 0) { qMax -= 1; } if (qMin < 3 || qMin > qMax) { continue; } if (qMin < lowerBig) { qMin = lowerBig | 1; } if (qMax >= upperExclusive) { qMax = upperExclusive - 1; if ((qMax & 1) == 0) { qMax -= 1; } } if (qMin > qMax || qMin > max64) { continue; } if (qMax > max64) { qMax = max64; if ((qMax & 1) == 0) { qMax -= 1; } } BigInteger stepBig = step64; BigInteger kMinBig = (qMin - 1 + (stepBig - 1)) / stepBig; BigInteger kMaxBig = (qMax - 1) / stepBig; BigInteger kStartBig = kMinBig - BigInteger.One; if (kStartBig < BigInteger.One) { kStartBig = kMinBig; } if (kStartBig < minK || kStartBig < effectiveResume) { kStartBig = BigInteger.Max(BigInteger.Max(minK, effectiveResume), kMinBig); } if (kStartBig > kMaxBig) { continue; } if (kMaxBig > ulong.MaxValue) { kMaxBig = ulong.MaxValue; } ulong kStart = (ulong)kStartBig; ulong kEnd = (ulong)kMaxBig; bool composite = CheckDivisors64Range( prime, step64, max64, kStart, kEnd, ref currentK, out bool processedRange, out foundDivisor); if (composite) { RecordPow2Hit(phase); FlushPow2State(currentK, phase); _currentPow2Phase = Pow2Phase.None; return true; } if (!processedRange) { processedAll = false; } } } if (currentK >= (ulong)effectiveResume) { FlushPow2State(currentK, phase); } return false; } finally { if (currentK >= (ulong)effectiveResume) { FlushPow2State(currentK, phase); } _currentPow2Phase = Pow2Phase.None; } } [MethodImpl(MethodImplOptions.AggressiveInlining)] private void EnsurePreparedSpecials(int specialRangePromille) { if (_preparedSpecialsInitialized && specialRangePromille == _preparedSpecialRangePromille) { return; } _preparedSpecialRangePromille = specialRangePromille; if (specialRangePromille == 0) { _preparedSpecialSingles = SpecialPercentPromille; _preparedSpecialRanges = Array.Empty<(int Start, int End)>(); _preparedSpecialsInitialized = true; return; } int[] specials = SpecialPercentPromille; if (specials.Length == 0) { _preparedSpecialSingles = Array.Empty(); _preparedSpecialRanges = Array.Empty<(int Start, int End)>(); _preparedSpecialsInitialized = true; return; } int[] sorted = new int[specials.Length]; Array.Copy(specials, sorted, specials.Length); Array.Sort(sorted); int range = specialRangePromille; List<(int Start, int End)> merged = new(sorted.Length); int currentStart = Math.Max(0, sorted[0] - range); int currentEnd = Math.Min(PercentScale, sorted[0] + range); for (int i = 1; i < sorted.Length; i++) { int start = Math.Max(0, sorted[i] - range); int end = Math.Min(PercentScale, sorted[i] + range); if (start <= currentEnd) { if (end > currentEnd) { currentEnd = end; } } else { merged.Add((currentStart, currentEnd)); currentStart = start; currentEnd = end; } } merged.Add((currentStart, currentEnd)); _preparedSpecialRanges = merged.ToArray(); _preparedSpecialSingles = Array.Empty(); _preparedSpecialsInitialized = true; } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private void FlushPow2State(BigInteger k, Pow2Phase phase) { if (phase == Pow2Phase.None || k < MinK) { return; } string targetPath = phase == Pow2Phase.Special ? _specialStateFilePath : _groupsStateFilePath; ref int counter = ref (phase == Pow2Phase.Special ? ref _specialStateCounter : ref _groupsStateCounter); ref BigInteger last = ref (phase == Pow2Phase.Special ? ref _specialLastSavedK : ref _groupsLastSavedK); if (k <= last) { return; } string? directory = Path.GetDirectoryName(targetPath); if (!string.IsNullOrEmpty(directory)) { Directory.CreateDirectory(directory); } File.AppendAllText(targetPath, k.ToString(CultureInfo.InvariantCulture) + Environment.NewLine); counter = 0; last = k; } #endif [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisors128Range( ulong prime, ushort primeDecimalMask, UInt128 step, BigInteger allowedMax, UInt128 startK, UInt128 endK, out BigInteger foundDivisor, out bool processedAll) { foundDivisor = BigInteger.Zero; processedAll = true; if (startK > endK) { return false; } UInt128 allowedMax128 = allowedMax <= BigIntegerNumbers.UInt128MaxValue ? (UInt128)allowedMax : UInt128.MaxValue; UInt128 maxKAllowed = (allowedMax128 - UInt128.One) / step; if (endK > maxKAllowed) { endK = maxKAllowed; if (startK > endK) { return false; } } UInt128 divisor = (step * startK) + UInt128.One; processedAll = true; for (UInt128 k = startK; k <= endK; k++) { ushort rem3 = (ushort)(divisor % 3); ushort rem5 = (ushort)(divisor % 5); ushort rem7 = (ushort)(divisor % 7); ushort rem11 = (ushort)(divisor % 11); ushort rem13 = (ushort)(divisor % 13); ushort rem17 = (ushort)(divisor % 17); ushort rem19 = (ushort)(divisor % 19); if (rem3 != 0 && rem5 != 0 && rem7 != 0 && rem11 != 0 && rem13 != 0 && rem17 != 0 && rem19 != 0) { int rem10 = (int)(divisor % 10); ushort rem8 = (ushort)(divisor % 8); if ((rem8 == 1 || rem8 == 7) && ((primeDecimalMask >> rem10) & 1) != 0) { BigInteger divisorBig = (BigInteger)divisor; RecordState((BigInteger)k); BigInteger powResult = BigInteger.ModPow(2, prime, divisorBig); if (powResult.IsOne) { foundDivisor = divisorBig; return true; } } } if (k == endK) { break; } divisor += step; } return false; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisorsLarge( ulong prime, ushort primeDecimalMask, BigInteger allowedMax, BigInteger currentK, BigInteger step, out bool processedAll, out BigInteger foundDivisor) { foundDivisor = BigInteger.Zero; if (allowedMax.IsZero) { processedAll = false; return false; } processedAll = true; BigInteger divisor = (step * currentK) + BigInteger.One; if (divisor > allowedMax) { return false; } BigInteger maxKAllowed = (allowedMax - BigInteger.One) / step; if (currentK > maxKAllowed) { return false; } BigInteger remainingIterationsBig = maxKAllowed - currentK + BigInteger.One; if (step <= ulong.MaxValue) { UInt128 step128 = (UInt128)(ulong)step; UInt128 maxKByStep = UInt128Numbers.MaxValueMinusOne / step128; bool allowedFits128 = allowedMax <= BigIntegerNumbers.UInt128MaxValue; if (allowedFits128) { UInt128 limitK = (((UInt128)allowedMax) - UInt128.One) / step128; if (limitK < maxKByStep) { maxKByStep = limitK; } } BigInteger maxKByStepBig = (BigInteger)maxKByStep; if (maxKByStepBig >= currentK) { if (currentK > BigIntegerNumbers.UInt128MaxValue) { currentK = maxKByStepBig + BigInteger.One; } else { UInt128 startK = (UInt128)currentK; bool composite128 = CheckDivisors128Range( prime, primeDecimalMask, step128, allowedMax, startK, maxKByStep, out BigInteger found128, out bool processed128); if (composite128) { foundDivisor = found128; return true; } if (!processed128) { processedAll = false; } currentK = maxKByStepBig + BigInteger.One; divisor = (step * currentK) + BigInteger.One; if (divisor > allowedMax) { return false; } } } } var rem10 = new BigIntegerCycleRemainderStepper(10); var rem8 = new BigIntegerCycleRemainderStepper(8); 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 = rem10.Initialize(divisor); BigInteger primeBig = (BigInteger)prime; while (remainingIterationsBig.Sign > 0) { bool passesSmallModuli = remainder10 == 1 || remainder10 == 7; if (passesSmallModuli && ((primeDecimalMask >> remainder10) & 1) != 0) { if (rem3.NextIsNotDivisible(divisor) && rem7.NextIsNotDivisible(divisor) && rem11.NextIsNotDivisible(divisor) && rem13.NextIsNotDivisible(divisor) && rem17.NextIsNotDivisible(divisor) && rem19.NextIsNotDivisible(divisor)) { byte remainder8 = rem8.ComputeNext(divisor); if (remainder8 == 1 || remainder8 == 7) { RecordState(currentK); BigInteger powResult = BigInteger.ModPow(BigIntegerNumbers.Two, primeBig, divisor); if (powResult.IsOne) { foundDivisor = divisor; processedAll = true; return true; } } } } currentK++; divisor += step; remainder10 = rem10.ComputeNext(divisor); remainingIterationsBig--; if (remainingIterationsBig.IsZero) { break; } } processedAll = true; foundDivisor = BigInteger.Zero; return false; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] 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; } ulong maxKAllowed = (allowedMax - 1UL) / step; if (endK > maxKAllowed) { endK = maxKAllowed; if (startK > endK) { 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); } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisors64( ulong prime, ulong step, ulong limit, ulong divisor, ref MersenneDivisorResidueStepper residueStepper, ref ulong currentK, out bool processedAll, out ulong foundDivisor) { #if DEVICE_GPU || DEVICE_HYBRID PrimeOrderCalculatorAccelerator gpu = Accelerator; #endif processedAll = true; foundDivisor = 0UL; ulong computedCycle; bool primeOrderFailed; ulong divisorCycle; bool computed; if (step > limit) { processedAll = true; foundDivisor = 0UL; return false; } if (divisor > limit) { processedAll = true; foundDivisor = 0UL; return false; } ulong maxKAllowed = ((limit - 1UL) / step); ulong targetK = ((divisor - 1UL) / step); if (targetK > maxKAllowed) { processedAll = true; foundDivisor = 0UL; return false; } ulong remainingIterations = (maxKAllowed - targetK) + 1UL; while (remainingIterations > 0UL) { if (residueStepper.IsAdmissible()) { MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); #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, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); #elif DEVICE_HYBRID divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( divisor, gpu, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); #else divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); #endif } else { divisorCycle = computedCycle; } RecordState(currentK); if (divisorCycle == prime) { foundDivisor = divisor; return true; } if (divisorCycle == 0UL) { Console.WriteLine($"Divisor cycle was not calculated for {prime}"); } } remainingIterations--; if (remainingIterations == 0UL) { return false; } divisor += step; currentK++; residueStepper.Advance(); } processedAll = true; foundDivisor = 0UL; return false; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private void RecordState(BigInteger k) { #if DivisorSet_Pow2Groups if (_currentPow2Phase == Pow2Phase.Special) { RecordStateInternal(k, _specialStateFilePath, ref _specialStateCounter, ref _specialLastSavedK, _currentPow2Phase); return; } if (_currentPow2Phase == Pow2Phase.Groups) { RecordStateInternal(k, _groupsStateFilePath, ref _groupsStateCounter, ref _groupsLastSavedK, _currentPow2Phase); return; } #endif RecordStateInternal(k, StateFilePath, ref _stateCounter, ref _lastSavedK, default); } #if DivisorSet_Pow2Groups private void RecordStateInternal(BigInteger k, string path, ref int counter, ref BigInteger lastSavedK, Pow2Phase phase) #else private void RecordStateInternal(BigInteger k, string path, ref int counter, ref BigInteger lastSavedK, object _) #endif { if (k <= lastSavedK) { return; } int next = counter + 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); counter = 0; lastSavedK = k; } else { counter = next; } } #if DivisorSet_Pow2Groups // Overload used from CheckDivisors64Range through RecordState with current phase routed automatically. private void RecordState(BigInteger k, Pow2Phase phase) { switch (phase) { case Pow2Phase.Special: RecordStateInternal(k, _specialStateFilePath, ref _specialStateCounter, ref _specialLastSavedK, phase); break; case Pow2Phase.Groups: RecordStateInternal(k, _groupsStateFilePath, ref _groupsStateCounter, ref _groupsLastSavedK, phase); break; default: RecordStateInternal(k, StateFilePath, ref _stateCounter, ref _lastSavedK, phase); break; } } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static void RecordPow2Hit(Pow2Phase phase) { long special = Interlocked.Read(ref _specialHits); long groups = Interlocked.Read(ref _groupHits); if (phase == Pow2Phase.Special) { special = Interlocked.Increment(ref _specialHits); } else if (phase == Pow2Phase.Groups) { groups = Interlocked.Increment(ref _groupHits); } Console.WriteLine($"By-divisor pow2groups hits: specials={special}, groups={groups}"); } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static int GetBitLength(BigInteger value) { if (value.IsZero) { return 0; } byte[] bytes = value.ToByteArray(); int lastIndex = bytes.Length - 1; return (lastIndex * 8) + BitOperations.Log2(bytes[lastIndex]) + 1; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool ProcessPow2PhaseBig( ulong prime, ushort primeDecimalMask, BigInteger allowedMax, BigInteger step, BigInteger resumeK, BigInteger minK, BigInteger segmentStartK, BigInteger segmentEndK, Pow2Phase phase, ReadOnlySpan<(int Start, int End)> ranges, ReadOnlySpan singles, ref bool processedAll, out BigInteger foundDivisor) { foundDivisor = BigInteger.Zero; if (allowedMax.IsZero || segmentStartK > segmentEndK) { return false; } int specialRangePromille = EnvironmentConfiguration.ByDivisorSpecialRange; if (specialRangePromille < 0) { specialRangePromille = 0; } specialRangePromille *= 100; // convert percent to promille-of-percentScale (10_000 = 100%) EnsurePreparedSpecials(specialRangePromille); BigInteger effectiveStart = segmentStartK; if (resumeK > effectiveStart) { effectiveStart = resumeK; } if (effectiveStart < minK) { effectiveStart = minK; } if (effectiveStart > segmentEndK) { return false; } BigInteger firstDivisor = (step * effectiveStart) + BigInteger.One; if (firstDivisor > allowedMax) { processedAll = false; return false; } int startBit = GetBitLength(firstDivisor) - 1; int maxBit = GetBitLength(allowedMax) - 1; _currentPow2Phase = phase; FlushPow2State(effectiveStart, phase); ReadOnlySpan effectiveSingles = singles; ReadOnlySpan<(int Start, int End)> effectiveRanges = ranges; if (phase == Pow2Phase.Special) { effectiveSingles = _preparedSpecialSingles; effectiveRanges = _preparedSpecialRanges; } try { for (int bit = startBit; bit <= maxBit; bit++) { BigInteger lower = BigInteger.One << bit; BigInteger upperExclusive = BigInteger.One << (bit + 1); if (lower > allowedMax) { break; } for (int i = 0; i < effectiveSingles.Length; i++) { int startProm = effectiveSingles[i]; int endProm = effectiveSingles[i]; if (specialRangePromille > 0) { startProm = Math.Max(0, effectiveSingles[i] - specialRangePromille); endProm = Math.Min(PercentScale, effectiveSingles[i] + specialRangePromille); } BigInteger scaledStart = (PercentScale + startProm) * lower; BigInteger scaledEnd = (PercentScale + (specialRangePromille == 0 ? Math.Min(PercentScale, startProm + 1) : endProm)) * lower; BigInteger qMin = (scaledStart + (PercentScale - 1)) / PercentScale; BigInteger qMax = scaledEnd / PercentScale; if (qMin.IsEven) { qMin += BigInteger.One; } if (qMax.IsEven) { qMax -= BigInteger.One; } if (qMin < 3 || qMin > qMax) { continue; } if (qMin < lower) { qMin = lower | BigInteger.One; } if (qMax >= upperExclusive) { qMax = upperExclusive - BigInteger.One; if ((qMax & BigInteger.One) == BigInteger.Zero) { qMax -= BigInteger.One; } } if (qMin > qMax || qMin > allowedMax) { continue; } if (qMax > allowedMax) { qMax = allowedMax; if ((qMax & BigInteger.One) == BigInteger.Zero) { qMax -= BigInteger.One; } } BigInteger kMin = (qMin - BigInteger.One + (step - BigInteger.One)) / step; BigInteger kMax = (qMax - BigInteger.One) / step; if (kMin < minK || kMin < effectiveStart) { kMin = BigInteger.Max(BigInteger.Max(minK, effectiveStart), kMin); } if (kMax > segmentEndK) { kMax = segmentEndK; } if (kMin > kMax) { continue; } BigInteger currentK = kMin; bool composite = CheckDivisorsBigRangePow2( prime, primeDecimalMask, step, allowedMax, kMin, kMax, phase, ref currentK, out bool processedRange, out BigInteger foundBig); if (composite) { RecordPow2Hit(phase); FlushPow2State(currentK, phase); _currentPow2Phase = Pow2Phase.None; foundDivisor = foundBig; return true; } if (!processedRange) { processedAll = false; } } for (int i = 0; i < effectiveRanges.Length; i++) { int startProm = effectiveRanges[i].Start; int endProm = effectiveRanges[i].End; if (specialRangePromille > 0 && phase == Pow2Phase.Special) { startProm = Math.Max(0, startProm - specialRangePromille); endProm = Math.Min(PercentScale, endProm + specialRangePromille); } BigInteger scaledStart = (PercentScale + startProm) * lower; BigInteger scaledEnd = (PercentScale + endProm) * lower; BigInteger qMin = (scaledStart + (PercentScale - 1)) / PercentScale; BigInteger qMax = scaledEnd / PercentScale; if ((qMin & BigInteger.One) == BigInteger.Zero) { qMin += BigInteger.One; } if ((qMax & BigInteger.One) == BigInteger.Zero) { qMax -= BigInteger.One; } if (qMin < 3 || qMin > qMax) { continue; } if (qMin < lower) { qMin = lower | BigInteger.One; } if (qMax >= upperExclusive) { qMax = upperExclusive - BigInteger.One; if ((qMax & BigInteger.One) == BigInteger.Zero) { qMax -= BigInteger.One; } } if (qMin > qMax || qMin > allowedMax) { continue; } if (qMax > allowedMax) { qMax = allowedMax; if ((qMax & BigInteger.One) == BigInteger.Zero) { qMax -= BigInteger.One; } } BigInteger kMin = (qMin - BigInteger.One + (step - BigInteger.One)) / step; BigInteger kMax = (qMax - BigInteger.One) / step; if (kMin < minK || kMin < effectiveStart) { kMin = BigInteger.Max(BigInteger.Max(minK, effectiveStart), kMin); } if (kMax > segmentEndK) { kMax = segmentEndK; } if (kMin > kMax) { continue; } BigInteger currentK = kMin; bool composite = CheckDivisorsBigRangePow2( prime, primeDecimalMask, step, allowedMax, kMin, kMax, phase, ref currentK, out bool processedRange, out BigInteger foundBig); if (composite) { RecordPow2Hit(phase); FlushPow2State(currentK, phase); _currentPow2Phase = Pow2Phase.None; foundDivisor = foundBig; return true; } if (!processedRange) { processedAll = false; } } } } finally { _currentPow2Phase = Pow2Phase.None; } return false; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisorsBigRangePow2( ulong prime, ushort primeDecimalMask, BigInteger step, BigInteger allowedMax, BigInteger startK, BigInteger endK, Pow2Phase phase, ref BigInteger currentK, out bool processedAll, out BigInteger foundDivisor) { if (startK < BigInteger.One || endK < startK) { processedAll = true; foundDivisor = BigInteger.Zero; return false; } currentK = startK; BigInteger divisor = (step * currentK) + BigInteger.One; if (divisor > allowedMax && !allowedMax.IsZero) { processedAll = false; foundDivisor = BigInteger.Zero; return false; } BigInteger iterations = (endK - currentK) + BigInteger.One; 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 (iterations > BigInteger.Zero) { bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((primeDecimalMask >> remainder10) & 1) != 0) { RecordState(currentK, phase); BigInteger powResult = BigInteger.ModPow(2, prime, divisor); if (powResult.IsOne) { foundDivisor = divisor; processedAll = true; return true; } } 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); iterations--; } processedAll = true; foundDivisor = BigInteger.Zero; return false; } #endif [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private void MarkPow2MinusOneChecked() { _pow2MinusOneChecked = true; if (!string.IsNullOrEmpty(_pow2MinusOneStateFilePath)) { try { string? directory = Path.GetDirectoryName(_pow2MinusOneStateFilePath); if (!string.IsNullOrEmpty(directory)) { Directory.CreateDirectory(directory); } File.WriteAllText(_pow2MinusOneStateFilePath, "checked"); } catch { // Best-effort persistence; ignore IO errors. } } } [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool TryFindPowerOfTwoMinusOneDivisor(ulong prime, BigInteger allowedMax, out BigInteger divisor) { divisor = BigInteger.Zero; if (_pow2MinusOneChecked) { return false; } if (!_checkedPow2MinusOne.TryAdd(prime, 0)) { return false; } BigInteger limit = allowedMax.IsZero ? new BigInteger(ulong.MaxValue) : allowedMax; if (limit <= ulong.MaxValue) { ulong limit64 = (ulong)limit; for (int x = 2; x < 64; x++) { ulong d = (1UL << x) - 1UL; if (d > limit64) { break; } if (BigInteger.ModPow(2, prime, d) == BigInteger.One) { divisor = d; return true; } } } BigInteger pow = BigIntegerNumbers.OneShiftedLeft63; for (int x = 63; ; x++) { pow <<= 1; BigInteger dBig = pow - BigInteger.One; if (dBig > limit) { break; } if (BigInteger.ModPow(2, prime, dBig) == BigInteger.One) { divisor = dBig; _checkedPow2MinusOne[prime] = 1; MarkPow2MinusOneChecked(); return true; } } _checkedPow2MinusOne[prime] = 1; MarkPow2MinusOneChecked(); return false; } }