using System.Numerics;
using System.Runtime.CompilerServices;
using System.Globalization;
using System.IO;
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<ulong> primes, Span<ulong> allowedMaxValues)
    {
        ulong divisorLimit = GetAllowedDivisorLimitForSpan(_divisorLimit);
        int length = primes.Length;
        for (int index = 0; index < length; index++)
        {
            allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit);
        }
    }

	public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu)
	{
		MersenneCpuDivisorScanSession? session = ThreadStaticPools.RentMersenneCpuDivisorSession();
		
		if (session is not null)
		{
			session.Configure(gpu, _orderDevice);
			return session;
		}

		return new MersenneCpuDivisorScanSession(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;

		// Keep the divisibility filters aligned with the divisor-cycle generator so the
		// CPU path never requests cycles that were skipped during cache creation.
		LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
		ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit);

		byte step10;
		byte step8;
		byte step3;
		byte step7;
		byte step11;
		byte step13;
		byte step17;
		byte step19;

		byte remainder10;
		byte remainder8;
		byte remainder3;
		byte remainder7;
		byte remainder11;
		byte remainder13;
		byte remainder17;
		byte remainder19;

		if (stepHigh == 0UL && limitHigh == 0UL)
		{
			step10 = (byte)(stepLow % 10UL);
			step8 = (byte)(stepLow & 7UL);
			step3 = (byte)(stepLow % 3UL);
			step7 = (byte)(stepLow % 7UL);
			step11 = (byte)(stepLow % 11UL);
			step13 = (byte)(stepLow % 13UL);
			step17 = (byte)(stepLow % 17UL);
			step19 = (byte)(stepLow % 19UL);

			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,
					decimalMask,
					step10,
					step8,
					step3,
					step7,
					step11,
					step13,
					step17,
					step19,
					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,
					decimalMask,
					step10,
					step8,
					step3,
					step7,
					step11,
					step13,
					step17,
					step19,
					ref currentK,
					out processedBottom,
					out foundDivisor);
			}

			processedAll = processedTop && processedBottom;
			return composite;
		}

		ulong divisorHigh = divisor.High;
		step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
		step8 = (byte)(stepLow % 8UL);
		step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL);
		step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL);
		step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL);
		step13 = (byte)((((stepHigh % 13UL) * 3UL) + (stepLow % 13UL)) % 13UL);
		step17 = (byte)(((stepHigh % 17UL) + (stepLow % 17UL)) % 17UL);
		step19 = (byte)((((stepHigh % 19UL) * 17UL) + (stepLow % 19UL)) % 19UL);

		remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL);
		remainder8 = (byte)(divisorLow % 8UL);
		remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL);
		remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
		remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL);
		remainder13 = (byte)((((divisorHigh % 13UL) * 3UL) + (divisorLow % 13UL)) % 13UL);
		remainder17 = (byte)(((divisorHigh % 17UL) + (divisorLow % 17UL)) % 17UL);
		remainder19 = (byte)((((divisorHigh % 19UL) * 17UL) + (divisorLow % 19UL)) % 19UL);

		var divisorPool = MontgomeryDivisorDataPool.Shared;
		while (divisor.CompareTo(limit) <= 0)
		{
			bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0;
			if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
			{
				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++;
			remainder10 = AddMod10(remainder10, step10);
			remainder8 = AddMod8(remainder8, step8);
			remainder3 = AddMod3(remainder3, step3);
			remainder7 = AddMod7(remainder7, step7);
			remainder11 = AddMod11(remainder11, step11);
			remainder13 = AddMod13(remainder13, step13);
			remainder17 = AddMod17(remainder17, step17);
			remainder19 = AddMod19(remainder19, step19);
		}

        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);

        byte step10 = (byte)(step % 10);
        byte step8 = (byte)(step % 8);
        byte step3 = (byte)(step % 3);
        byte step7 = (byte)(step % 7);
        byte step11 = (byte)(step % 11);
        byte step13 = (byte)(step % 13);
        byte step17 = (byte)(step % 17);
        byte step19 = (byte)(step % 19);

        byte remainder10 = (byte)(divisor % 10);
        byte remainder8 = (byte)(divisor % 8);
        byte remainder3 = (byte)(divisor % 3);
        byte remainder7 = (byte)(divisor % 7);
        byte remainder11 = (byte)(divisor % 11);
        byte remainder13 = (byte)(divisor % 13);
        byte remainder17 = (byte)(divisor % 17);
        byte remainder19 = (byte)(divisor % 19);

        while (divisor <= allowedMax)
        {
            bool passesSmallModuli = remainder3 != 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 = AddMod10(remainder10, step10);
            remainder8 = AddMod8(remainder8, step8);
            remainder3 = AddMod3(remainder3, step3);
            remainder7 = AddMod7(remainder7, step7);
            remainder11 = AddMod11(remainder11, step11);
            remainder13 = AddMod13(remainder13, step13);
            remainder17 = AddMod17(remainder17, step17);
            remainder19 = AddMod19(remainder19, step19);

            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,
		ushort decimalMask,
		byte step10,
		byte step8,
		byte step3,
		byte step7,
		byte step11,
		byte step13,
		byte step17,
		byte step19,
		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;

		byte remainder10 = (byte)(startDivisor % 10UL);
		byte remainder8 = (byte)(startDivisor & 7UL);
		byte remainder3 = (byte)(startDivisor % 3UL);
		byte remainder7 = (byte)(startDivisor % 7UL);
		byte remainder11 = (byte)(startDivisor % 11UL);
		byte remainder13 = (byte)(startDivisor % 13UL);
		byte remainder17 = (byte)(startDivisor % 17UL);
		byte remainder19 = (byte)(startDivisor % 19UL);

		return CheckDivisors64(
			gpu,
			prime,
			step,
			rangeLimit,
			startDivisor,
			decimalMask,
			step10,
			step8,
			step3,
			step7,
			step11,
			step13,
			step17,
			step19,
			remainder10,
			remainder8,
			remainder3,
			remainder7,
			remainder11,
			remainder13,
			remainder17,
			remainder19,
			ref currentK,
			out processedAll,
			out foundDivisor);
	}

private bool CheckDivisors64(
		PrimeOrderCalculatorAccelerator gpu,
		ulong prime,
		ulong step,
		ulong limit,
		ulong divisor,
		ushort decimalMask,
		byte step10,
		byte step8,
		byte step3,
		byte step7,
		byte step11,
		byte step13,
		byte step17,
		byte step19,
		byte remainder10,
		byte remainder8,
		byte remainder3,
		byte remainder7,
		byte remainder11,
		byte remainder13,
		byte remainder17,
		byte remainder19,
		ref ulong currentK,
		out bool processedAll,
		out ulong foundDivisor)
	{
		bool canAdvance = step <= limit;
		Queue<MontgomeryDivisorData> divisorPool = MontgomeryDivisorDataPool.Shared;
		if (!canAdvance)
		{
			if (divisor <= limit)
			{
				if (remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0 && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
				{
					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)
		{
			bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0;
			if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0)
			{
				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++;
			remainder10 = AddMod10(remainder10, step10);
			remainder8 = AddMod8(remainder8, step8);
			remainder3 = AddMod3(remainder3, step3);
			remainder7 = AddMod7(remainder7, step7);
			remainder11 = AddMod11(remainder11, step11);
			remainder13 = AddMod13(remainder13, step13);
			remainder17 = AddMod17(remainder17, step17);
			remainder19 = AddMod19(remainder19, step19);
		}
	}

	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<int> 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;
            }

            BigInteger x = BigInteger.ModPow(baseValue, d, value);
            if (x == 1 || x == value - 1)
            {
                continue;
            }

            bool witnessFound = true;
            for (int r = 1; r < s; r++)
            {
                x = BigInteger.ModPow(x, 2, value);
                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;
    }


}