// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// IntPrime.cpp
// Primality testing implementation

#include "math/core/mp/Int/IntPrime.hpp"
#include "math/core/mp/Int/IntOps.hpp"
#include "math/core/mp/Int/IntModular.hpp"
#include "math/core/mp/Int/IntGCD.hpp"
#include "math/core/mp/Int/IntSqrt.hpp"
#include "math/core/mp/ThreadPool.hpp"
#include <random>
#include <vector>
#include <numeric>
#include <cstdint>
#include <atomic>
#include <mutex>

namespace sangi {

// ============================================================================
// Building the table for single-limb batch GCD (executed only once, on first call)
// ============================================================================

const std::vector<IntPrime::PrimeBatch>& IntPrime::getBatches() {
    static std::vector<PrimeBatch> batches = []() {
        std::vector<PrimeBatch> b;
        PrimeBatch current;
        current.product = 1;
        current.start_index = 0;
        current.count = 0;

        for (int i = 0; i < SMALL_PRIMES_COUNT; i++) {
            uint64_t p = static_cast<uint64_t>(SMALL_PRIMES[i]);
            if (current.count > 0 && current.product > UINT64_MAX / p) {
                // Would overflow, so start a new batch
                b.push_back(current);
                current.product = p;
                current.start_index = i;
                current.count = 1;
            } else {
                current.product *= p;
                current.count++;
            }
        }
        if (current.count > 0) {
            b.push_back(current);
        }
        // Precompute mpn::divmod_1 preinv (avoids calling invert_limb every time)
        for (auto& batch : b) {
            batch.shift = static_cast<unsigned>(std::countl_zero(batch.product));
            batch.product_norm = batch.product << batch.shift;
            batch.product_dinv = mpn::invert_limb(batch.product_norm);
        }
        return b;
    }();
    return batches;
}

// ============================================================================
// isDivisibleBySmallPrime — single-limb batch GCD version
// Assumes n is odd (the caller has already checked for even values)
// ============================================================================

// MASM ASM function (Win MSVC + BMI2/ADX, mpn_x64_mod1.asm)
// Equivalent to GMP's mpn_mod_1_1p: computes 1-limb mod via Möller-Granlund preinv.
#if defined(SANGI_INT_HAS_ASM)
extern "C" uint64_t mpn_mod_1_preinv_asm(const uint64_t* ap, uint64_t an,
                                         uint64_t d_norm, uint64_t dinv,
                                         uint64_t shift);
#endif

// C++ fallback (for no-ASM builds)
static inline uint64_t mod_1_preinv_cxx(const uint64_t* a, size_t an,
                                        uint64_t d_norm, uint64_t dinv,
                                        unsigned shift) {
    if (an == 0) return 0;
    if (shift == 0) {
        uint64_t r = 0;
        for (size_t i = an; i-- > 0; ) {
            auto [qq, rr] = mpn::udiv_qrnnd_preinv(r, a[i], d_norm, dinv);
            (void)qq;
            r = rr;
        }
        return r;
    }
    unsigned rshift = 64 - shift;
    uint64_t r = a[an - 1] >> rshift;
    for (size_t i = an; i-- > 1; ) {
        uint64_t n1 = (a[i] << shift) | (a[i - 1] >> rshift);
        auto [qq, rr] = mpn::udiv_qrnnd_preinv(r, n1, d_norm, dinv);
        (void)qq;
        r = rr;
    }
    {
        uint64_t n1 = a[0] << shift;
        auto [qq, rr] = mpn::udiv_qrnnd_preinv(r, n1, d_norm, dinv);
        (void)qq;
        r = rr;
    }
    return r >> shift;  // Denormalize the remainder
}

// preinv version that computes only the remainder (does not write q; shift/dinv are precomputed).
// Delegates to MASM ASM if available, otherwise uses the C++ fallback.
static inline uint64_t mod_1_preinv(const uint64_t* a, size_t an,
                                    uint64_t d_norm, uint64_t dinv,
                                    unsigned shift) {
#if defined(SANGI_INT_HAS_ASM)
    return mpn_mod_1_preinv_asm(a, static_cast<uint64_t>(an), d_norm, dinv,
                                static_cast<uint64_t>(shift));
#else
    return mod_1_preinv_cxx(a, an, d_norm, dinv, shift);
#endif
}

bool IntPrime::isDivisibleBySmallPrime(const Int& n) {
    const auto& batches = getBatches();

    // raw mpn path: avoids Int copy/heap alloc and computes 1-limb mod directly
    // via preinv mpn::mod_1 (the old Int operator% path takes ~292ns for 1024-bit n).
    // shift/d_norm/dinv are pre-computed in PrimeBatch (see getBatches).
    const uint64_t* nd = n.data();
    size_t nn = n.size();
    if (nn == 0) return false;

    for (const auto& batch : batches) {
        uint64_t r = mod_1_preinv(nd, nn, batch.product_norm, batch.product_dinv, batch.shift);
        if (r == 0) {
            return true;  // The entire batch.product divides n
        }
        uint64_t g = std::gcd(r, batch.product);
        if (g > 1) {
            // n itself might be a prime in the small-prime table (equality check for 1-limb)
            if (nn == 1) {
                for (int j = 0; j < batch.count; j++) {
                    if (nd[0] == static_cast<uint64_t>(SMALL_PRIMES[batch.start_index + j])) {
                        return false;
                    }
                }
            }
            return true;
        }
    }
    return false;
}

// ============================================================================
// extractSmallPrimeFactors — extract small-prime factors in bulk via single-limb batch GCD
// remaining: input is odd (the prime factor 2 has been removed beforehand), output is the value with small-prime factors removed
// factors: appends the factors found as {prime, exponent}
// ============================================================================

void IntPrime::extractSmallPrimeFactors(
    Int& remaining,
    std::vector<std::pair<Int, int>>& factors)
{
    if (remaining.isOne()) return;

    const auto& batches = getBatches();

    // Process the batches of larger primes first
    // Dividing by large prime factors first reduces the digit count of remaining sooner,
    // making N mod batch_product faster for subsequent batches
    for (int bi = static_cast<int>(batches.size()) - 1; bi >= 0; bi--) {
        const auto& batch = batches[bi];
        if (remaining.isOne()) break;

        // N mod batch.product (single-limb division)
        Int r = remaining % Int(batch.product);

        uint64_t g;
        if (r.isZero()) {
            g = batch.product;
        } else {
            uint64_t rv = r.toUInt64();
            g = std::gcd(rv, batch.product);
        }

        if (g <= 1) continue;  // None of the primes in this batch divide N

        // g > 1 → individually check which primes in this batch divide remaining
        // Dividing by larger prime factors first reduces the digit count of remaining sooner
        for (int j = batch.count - 1; j >= 0; j--) {
            if (remaining.isOne()) break;

            int p = SMALL_PRIMES[batch.start_index + j];
            if (g % p != 0) continue;  // This prime does not divide N

            // Divide out p as many times as possible
            Int prime(p);
            int exp = 0;
            while (true) {
                Int quotient = remaining / prime;
                Int rem = remaining % prime;
                if (!rem.isZero()) break;
                remaining = quotient;
                exp++;
            }

            if (exp > 0) {
                factors.push_back({Int(p), exp});
            }
        }
    }
}

// ============================================================================
// Miller-Rabin primality test
// ============================================================================

bool IntPrime::isMillerRabinPrime(const Int& n, int k) {
    // 1. Check special states
    if (n.isNaN() || n.isInfinite()) {
        return false;
    }

    // 2. Values less than 2 are not prime
    if (n < 2) {
        return false;
    }

    // 3. 2 is prime
    if (n == 2) {
        return true;
    }

    // 4. Even numbers are not prime (except 2)
    if (!n.getBit(0)) {
        return false;
    }

    // 5. If the value is within the small-prime table, decide directly
    //    The maximum value in SMALL_PRIMES is 1987
    if (n <= 1987) {
        for (int i = 0; i < SMALL_PRIMES_COUNT; i++) {
            if (n == Int(SMALL_PRIMES[i])) {
                return true;
            }
        }
        // An odd number ≤ 1987 not in the small-prime list → composite
        return false;
    }

    // 6. Divisibility check by small primes via batch GCD
    if (isDivisibleBySmallPrime(n)) {
        return false;
    }

    // 7. Decompose n-1 = 2^s × d (d is odd)
    Int n_minus_1 = n - 1;
    Int d = n_minus_1;
    int s = static_cast<int>(d.countTrailingZeros());
    d = d >> s;

    // 8. k rounds of the Miller-Rabin test
    std::random_device rd;
    std::mt19937_64 gen(rd());

    // Use small primes as bases, and random ones if there are not enough
    for (int round = 0; round < k; round++) {
        Int a;
        if (round < SMALL_PRIMES_COUNT) {
            a = Int(SMALL_PRIMES[round]);
        } else {
            a = Int(2 + static_cast<int64_t>(gen() % 20));
        }

        if (a >= n) {
            a = IntModular::mod(a, n);
        }
        // Skip degenerate witnesses (0 or 1)
        if (a <= Int::One()) {
            continue;
        }

        Int x = IntModular::powerMod(a, d, n);

        if (x.isOne()) {
            continue;
        }

        if (x == n_minus_1) {
            continue;
        }

        bool found_minus_one = false;
        for (int r = 1; r < s; r++) {
            x = IntModular::powerMod(x, 2, n);

            if (x == n_minus_1) {
                found_minus_one = true;
                break;
            }
        }

        if (!found_minus_one) {
            return false;
        }
    }

    return true;
}

// ============================================================================
// Fermat primality test
// ============================================================================

bool IntPrime::isFermatPrime(const Int& n, int k) {
    if (n.isNaN() || n.isInfinite()) {
        return false;
    }

    if (n < 2) {
        return false;
    }

    if (n == 2) {
        return true;
    }

    if (!n.getBit(0)) {
        return false;
    }

    if (n <= 1987) {
        for (int i = 0; i < SMALL_PRIMES_COUNT; i++) {
            if (n == Int(SMALL_PRIMES[i])) {
                return true;
            }
        }
        return false;
    }

    if (isDivisibleBySmallPrime(n)) {
        return false;
    }

    Int n_minus_1 = n - 1;

    for (int i = 0; i < k && i < SMALL_PRIMES_COUNT; i++) {
        Int a(SMALL_PRIMES[i]);

        if (!IntGCD::gcd(a, n).isOne()) {
            continue;
        }

        Int result = IntModular::powerMod(a, n_minus_1, n);

        if (!result.isOne()) {
            return false;
        }
    }

    return true;
}

// ============================================================================
// nextPrime
// ============================================================================

Int IntPrime::nextPrime(const Int& n, int k) {
    if (n.isNaN()) {
        return Int::NaN();
    }
    if (n.isInfinite()) {
        if (n.getSign() > 0) {
            return Int::PositiveInfinity();
        } else {
            return Int(2);
        }
    }

    if (n < 2) {
        return Int(2);
    }

    Int candidate = n + 1;

    if (candidate > 2 && candidate.isEven()) {
        candidate += 1;
    }

    while (true) {
        if (isMillerRabinPrime(candidate, k)) {
            return candidate;
        }
        if (candidate == 2) {
            candidate = Int(3);
        } else {
            candidate += 2;
        }
    }
}

// ============================================================================
// prevPrime
// ============================================================================

Int IntPrime::prevPrime(const Int& n, int k) {
    if (n.isNaN()) {
        return Int::NaN();
    }
    if (n.isInfinite()) {
        return Int::NaN();
    }

    if (n <= 2) {
        return Int::NaN();
    }
    if (n <= 3) {
        return Int(2);
    }

    Int candidate = n - 1;

    if (candidate.isEven()) {
        candidate -= 1;
    }

    while (candidate >= 2) {
        if (isMillerRabinPrime(candidate, k)) {
            return candidate;
        }
        candidate -= 2;

        if (candidate < 3) {
            candidate = Int(2);
        }
    }

    return Int::NaN();
}

// ============================================================================
// findFactor_Fermat
// ============================================================================

Int IntPrime::findFactor_Fermat(const Int& n, uint64_t max_iterations) {
    // If n is even, return 2
    if (n.isEven()) {
        return Int(2);
    }
    if (n <= 1) {
        return n;
    }

    // a = ceil(sqrt(n))
    Int a = IntOps::sqrt(n);
    if (a * a == n) {
        return a;  // Perfect square
    }
    a += 1;

    // Compute b² = a² - n, and increment a until b² is a perfect square
    Int b2 = a * a - n;

    for (uint64_t iter = 0; iter < max_iterations; ++iter) {
        Int b = IntOps::sqrt(b2);
        if (b * b == b2) {
            // n = (a+b)(a-b)
            Int factor = a - b;
            if (factor > 1 && factor < n) {
                return factor;
            }
            factor = a + b;
            if (factor > 1 && factor < n) {
                return factor;
            }
        }
        // a++ → b2 = (a+1)² - n = a² - n + 2a + 1 = b2 + 2a + 1
        b2 = b2 + a + a + 1;
        a += 1;
    }

    return n;  // No factor found
}

// ============================================================================
// Pollard rho — fast uint64-only version (n ≤ 63 bits)
// Eliminates Int construction entirely; sped up with 128-bit multiply + binary GCD
// ============================================================================

namespace {

// 64-bit modular multiply: (a * b) mod m (a, b < m < 2^63)
inline uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t m) {
#ifdef _MSC_VER
    uint64_t hi;
    uint64_t lo = _umul128(a, b, &hi);
    uint64_t rem;
    _udiv128(hi, lo, m, &rem);
    return rem;
#else
    return static_cast<uint64_t>(
        (static_cast<unsigned __int128>(a) * b) % m);
#endif
}

// 64-bit binary GCD
inline uint64_t gcd64(uint64_t a, uint64_t b) {
    if (a == 0) return b;
    if (b == 0) return a;
#ifdef _MSC_VER
    unsigned long shift;
    _BitScanForward64(&shift, a | b);
    unsigned long sa;
    _BitScanForward64(&sa, a);
    a >>= sa;
#else
    int shift = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);
#endif
    do {
#ifdef _MSC_VER
        unsigned long sb;
        _BitScanForward64(&sb, b);
        b >>= sb;
#else
        b >>= __builtin_ctzll(b);
#endif
        if (a > b) { uint64_t t = a; a = b; b = t; }
        b -= a;
    } while (b != 0);
    return a << shift;
}

// Brent's rho — uint64-only, batch GCD (m=128)
uint64_t pollardRho64(uint64_t n) {
    if (n % 2 == 0) return 2;
    if (n < 4) return n;

    for (uint64_t c = 1; c <= 20; ++c) {
        uint64_t x = 2, y = 2, d = 1;
        uint64_t ys = 0, q = 1;
        constexpr uint64_t batch_size = 128;
        uint64_t r = 1;

        while (d == 1) {
            x = y;
            for (uint64_t i = 0; i < r; ++i)
                y = mulmod64(y, y, n) + c;  // c < n, so no overflow (n < 2^63)

            uint64_t k = 0;
            while (k < r && d == 1) {
                ys = y;
                uint64_t b = (r - k < batch_size) ? (r - k) : batch_size;
                for (uint64_t i = 0; i < b; ++i) {
                    y = mulmod64(y, y, n) + c;
                    uint64_t diff = (y > x) ? (y - x) : (x - y);
                    q = mulmod64(q, diff, n);
                }
                d = gcd64(q, n);
                k += b;
            }
            r *= 2;
            if (r > 1000000) break;
        }

        if (d == n) {
            // Backtrack
            d = 1;
            while (d == 1) {
                ys = mulmod64(ys, ys, n) + c;
                uint64_t diff = (ys > x) ? (ys - x) : (x - ys);
                d = gcd64(diff, n);
            }
        }

        if (d > 1 && d < n) return d;
    }
    return n;
}

// ============================================================================
// Pollard rho — 128-bit Montgomery-only version (64 < n ≤ 127 bits)
// ============================================================================

struct Mont128 {
    uint64_t n_lo, n_hi;     // n (128-bit)
    uint64_t inv_lo;         // -n^{-1} mod 2^64 (low part of the Montgomery inverse)

    // 128-bit addition: a + b, subtract n if the result is ≥ n
    static void add128(uint64_t& r_lo, uint64_t& r_hi,
                       uint64_t a_lo, uint64_t a_hi,
                       uint64_t b_lo, uint64_t b_hi) {
        r_lo = a_lo + b_lo;
        r_hi = a_hi + b_hi + (r_lo < a_lo ? 1 : 0);
    }

    // 128-bit subtraction: a - b (assumes a >= b)
    static void sub128(uint64_t& r_lo, uint64_t& r_hi,
                       uint64_t a_lo, uint64_t a_hi,
                       uint64_t b_lo, uint64_t b_hi) {
        r_hi = a_hi - b_hi - (a_lo < b_lo ? 1 : 0);
        r_lo = a_lo - b_lo;
    }

    // 128-bit comparison: a >= b
    static bool ge128(uint64_t a_lo, uint64_t a_hi,
                      uint64_t b_lo, uint64_t b_hi) {
        return (a_hi > b_hi) || (a_hi == b_hi && a_lo >= b_lo);
    }

    // a == 0
    static bool isZero128(uint64_t lo, uint64_t hi) {
        return lo == 0 && hi == 0;
    }

    // Montgomery REDC: T (256-bit) → T * R^{-1} mod n  (R = 2^128)
    // 2-limb SOS method: add m_i * n to T to cancel the low limbs
    void redc(uint64_t& r_lo, uint64_t& r_hi,
              uint64_t t0, uint64_t t1, uint64_t t2, uint64_t t3) const {
        // === Step 1: m0 = t0 * inv_lo mod 2^64, T += m0 * n ===
        uint64_t m0 = t0 * inv_lo;
        // m0 * n = m0*n_lo + m0*n_hi * 2^64  (up to 192 bits)
        uint64_t p0_hi, p0_lo = _umul128(m0, n_lo, &p0_hi);
        uint64_t p1_hi, p1_lo = _umul128(m0, n_hi, &p1_hi);

        // T[0] += p0_lo → 0 by construction, carry = (the high bit of T[0] + p0_lo)
        uint64_t c = 0;
        {
            uint64_t s = t0 + p0_lo;
            c = (s < t0) ? 1ULL : 0ULL;
            // s == 0 by construction (m0 * n_lo ≡ -t0 mod 2^64)
        }
        // T[1] += p0_hi + p1_lo + c
        {
            uint64_t s = t1 + p0_hi;
            uint64_t c1 = (s < t1) ? 1ULL : 0ULL;
            uint64_t s2 = s + p1_lo;
            c1 += (s2 < s) ? 1ULL : 0ULL;
            uint64_t s3 = s2 + c;
            c1 += (s3 < s2) ? 1ULL : 0ULL;
            t1 = s3;  // The new t1 (cancelled in the next step)
            c = c1;
        }
        // T[2] += p1_hi + c
        {
            uint64_t s = t2 + p1_hi;
            uint64_t c1 = (s < t2) ? 1ULL : 0ULL;
            uint64_t s2 = s + c;
            c1 += (s2 < s) ? 1ULL : 0ULL;
            t2 = s2;
            c = c1;
        }
        t3 += c;

        // === Step 2: m1 = t1 * inv_lo mod 2^64, T += m1 * n * 2^64 ===
        uint64_t m1 = t1 * inv_lo;
        uint64_t q0_hi, q0_lo = _umul128(m1, n_lo, &q0_hi);
        uint64_t q1_hi, q1_lo = _umul128(m1, n_hi, &q1_hi);

        c = 0;
        // T[1] += q0_lo → 0 by construction
        {
            uint64_t s = t1 + q0_lo;
            c = (s < t1) ? 1ULL : 0ULL;
        }
        // T[2] += q0_hi + q1_lo + c
        {
            uint64_t s = t2 + q0_hi;
            uint64_t c1 = (s < t2) ? 1ULL : 0ULL;
            uint64_t s2 = s + q1_lo;
            c1 += (s2 < s) ? 1ULL : 0ULL;
            uint64_t s3 = s2 + c;
            c1 += (s3 < s2) ? 1ULL : 0ULL;
            r_lo = s3;
            c = c1;
        }
        // T[3] += q1_hi + c
        {
            uint64_t s = t3 + q1_hi;
            uint64_t c1 = (s < t3) ? 1ULL : 0ULL;
            uint64_t s2 = s + c;
            c1 += (s2 < s) ? 1ULL : 0ULL;
            r_hi = s2;
            c = c1;
        }

        // Conditional subtraction: if result >= n then result -= n
        if (c > 0 || ge128(r_lo, r_hi, n_lo, n_hi)) {
            sub128(r_lo, r_hi, r_lo, r_hi, n_lo, n_hi);
        }
    }

    // Montgomery multiplication: a * b * R^{-1} mod n
    void mont_mul(uint64_t& r_lo, uint64_t& r_hi,
                  uint64_t a_lo, uint64_t a_hi,
                  uint64_t b_lo, uint64_t b_hi) const {
        // a * b → 256-bit (t3:t2:t1:t0)
        // = a_lo*b_lo + (a_lo*b_hi + a_hi*b_lo)<<64 + a_hi*b_hi<<128
        uint64_t p0_hi, p0_lo = _umul128(a_lo, b_lo, &p0_hi);
        uint64_t p1_hi, p1_lo = _umul128(a_lo, b_hi, &p1_hi);
        uint64_t p2_hi, p2_lo = _umul128(a_hi, b_lo, &p2_hi);
        uint64_t p3_hi, p3_lo = _umul128(a_hi, b_hi, &p3_hi);

        uint64_t t0 = p0_lo;
        uint64_t t1 = p0_hi;
        uint64_t t2 = p3_lo;
        uint64_t t3 = p3_hi;

        // t1 += p1_lo + p2_lo
        uint64_t c = 0;
        t1 += p1_lo; c += (t1 < p1_lo) ? 1 : 0;
        t1 += p2_lo; c += (t1 < p2_lo) ? 1 : 0;
        // t2 += p1_hi + p2_hi + c
        t2 += p1_hi; uint64_t c2 = (t2 < p1_hi) ? 1 : 0;
        t2 += p2_hi; c2 += (t2 < p2_hi) ? 1 : 0;
        t2 += c;     c2 += (t2 < c) ? 1 : 0;
        t3 += c2;

        redc(r_lo, r_hi, t0, t1, t2, t3);
    }
};

// 128-bit binary GCD → uint64_t result (assumes the GCD fits in 64 bits)
uint64_t gcd128_to64(uint64_t a_lo, uint64_t a_hi,
                     uint64_t b_lo, uint64_t b_hi) {
    // Simple version: use Int's GCD
    // In rho's batch GCD, this is the GCD of n and the 128-bit value of product reduced mod n
    // The result is a factor, so it is at most 64 bits
    if (a_hi == 0 && a_lo == 0) {
        if (b_hi == 0) return b_lo;
        // b is large → use Int
    }
    if (b_hi == 0 && b_lo == 0) {
        if (a_hi == 0) return a_lo;
    }
    // Fast path if both are 64-bit
    if (a_hi == 0 && b_hi == 0) return gcd64(a_lo, b_lo);

    // Fallback via Int
    Int ai, bi;
    if (a_hi == 0) ai = Int(a_lo);
    else { ai = Int(a_hi); ai <<= 64; ai += a_lo; }
    if (b_hi == 0) bi = Int(b_lo);
    else { bi = Int(b_hi); bi <<= 64; bi += b_lo; }
    Int g = IntGCD::gcd(ai, bi);
    return g.toUInt64();
}

// Brent's rho — 128-bit Montgomery version
// n is 2 limbs (64 < bitLength <= 127, odd)
uint64_t pollardRho128(uint64_t n_lo, uint64_t n_hi) {
    // Montgomery setup
    Mont128 mont;
    mont.n_lo = n_lo;
    mont.n_hi = n_hi;

    // -n^{-1} mod 2^64 (Newton's method: x_{i+1} = x_i * (2 - n * x_i))
    // Initial value x_0 = 1 (for odd n, n*1 ≡ 1 mod 2)
    // Converges to 64-bit precision in 6 iterations
    uint64_t inv = 1;
    for (int i = 0; i < 6; ++i)
        inv = inv * (2 - n_lo * inv);
    // inv = n^{-1} mod 2^64; negate to get -n^{-1} mod 2^64
    mont.inv_lo = ~inv + 1;  // = -inv mod 2^64

    // Compute R = 2^128 mod n, R² mod n (initialized with ordinary division)
    // Simple: compute R² mod n via Int
    Int nInt;
    nInt = Int(n_hi); nInt <<= 64; nInt += n_lo;
    Int R = Int(1); R <<= 128; R = R % nInt;
    Int R2 = R * R % nInt;

    uint64_t R_lo, R_hi, R2_lo, R2_hi;
    if (R.bitLength() <= 64) { R_lo = R.toUInt64(); R_hi = 0; }
    else {
        Int Rhi = R >> 64;
        R_lo = (R - (Rhi << 64)).toUInt64();
        R_hi = Rhi.toUInt64();
    }
    if (R2.bitLength() <= 64) { R2_lo = R2.toUInt64(); R2_hi = 0; }
    else {
        Int R2hi = R2 >> 64;
        R2_lo = (R2 - (R2hi << 64)).toUInt64();
        R2_hi = R2hi.toUInt64();
    }

    // toMont(x) = x * R2 (Montgomery multiplication)
    // fromMont(x) = x * 1 (REDC)

    for (uint64_t c_val = 1; c_val <= 20; ++c_val) {
        // Convert c, x, y to Montgomery form
        uint64_t c_lo, c_hi;
        mont.mont_mul(c_lo, c_hi, c_val, 0, R2_lo, R2_hi);  // c * R mod n
        uint64_t two_lo, two_hi;
        mont.mont_mul(two_lo, two_hi, 2, 0, R2_lo, R2_hi);    // 2 * R mod n

        uint64_t x_lo = two_lo, x_hi = two_hi;
        uint64_t y_lo = two_lo, y_hi = two_hi;
        uint64_t d = 1;
        uint64_t ys_lo = 0, ys_hi = 0;
        // q = 1 in Montgomery form = R mod n
        uint64_t q_lo = R_lo, q_hi = R_hi;
        constexpr uint64_t batch_size = 128;
        uint64_t r = 1;

        while (d == 1) {
            x_lo = y_lo; x_hi = y_hi;
            for (uint64_t i = 0; i < r; ++i) {
                // y = y² + c (mod n) in Montgomery form
                uint64_t t_lo, t_hi;
                mont.mont_mul(t_lo, t_hi, y_lo, y_hi, y_lo, y_hi);
                Mont128::add128(y_lo, y_hi, t_lo, t_hi, c_lo, c_hi);
                if (Mont128::ge128(y_lo, y_hi, n_lo, n_hi))
                    Mont128::sub128(y_lo, y_hi, y_lo, y_hi, n_lo, n_hi);
            }

            uint64_t k = 0;
            while (k < r && d == 1) {
                ys_lo = y_lo; ys_hi = y_hi;
                uint64_t b = (r - k < batch_size) ? (r - k) : batch_size;
                for (uint64_t i = 0; i < b; ++i) {
                    // y = y² + c
                    uint64_t t_lo, t_hi;
                    mont.mont_mul(t_lo, t_hi, y_lo, y_hi, y_lo, y_hi);
                    Mont128::add128(y_lo, y_hi, t_lo, t_hi, c_lo, c_hi);
                    if (Mont128::ge128(y_lo, y_hi, n_lo, n_hi))
                        Mont128::sub128(y_lo, y_hi, y_lo, y_hi, n_lo, n_hi);

                    // diff = |x - y| (take the difference while still in Montgomery form)
                    uint64_t diff_lo, diff_hi;
                    if (Mont128::ge128(y_lo, y_hi, x_lo, x_hi))
                        Mont128::sub128(diff_lo, diff_hi, y_lo, y_hi, x_lo, x_hi);
                    else
                        Mont128::sub128(diff_lo, diff_hi, x_lo, x_hi, y_lo, y_hi);

                    // q = q * diff (Montgomery multiplication)
                    uint64_t nq_lo, nq_hi;
                    mont.mont_mul(nq_lo, nq_hi, q_lo, q_hi, diff_lo, diff_hi);
                    q_lo = nq_lo; q_hi = nq_hi;
                }
                // GCD(q, n) — convert q back to normal form, then take GCD
                uint64_t qq_lo, qq_hi;
                mont.redc(qq_lo, qq_hi, q_lo, q_hi, 0, 0);  // fromMont(q)
                d = gcd128_to64(qq_lo, qq_hi, n_lo, n_hi);
                k += b;
            }
            r *= 2;
            if (r > 1000000) break;
        }

        if (d == n_lo && n_hi == 0) {
            // Case d == n: backtrack (a simple comparison is not possible when n_hi != 0)
            goto backtrack;
        }
        if (n_hi != 0) {
            // d is 64-bit but n is 128-bit → d == n cannot happen
            // However, if d equals n_lo and n_hi is nonzero, then d < n is guaranteed
        }
        // Check for d == n (128-bit)
        if (d == n_lo && n_hi == 0) {
            backtrack:
            d = 1;
            while (d == 1) {
                uint64_t t_lo, t_hi;
                mont.mont_mul(t_lo, t_hi, ys_lo, ys_hi, ys_lo, ys_hi);
                Mont128::add128(ys_lo, ys_hi, t_lo, t_hi, c_lo, c_hi);
                if (Mont128::ge128(ys_lo, ys_hi, n_lo, n_hi))
                    Mont128::sub128(ys_lo, ys_hi, ys_lo, ys_hi, n_lo, n_hi);

                uint64_t diff_lo, diff_hi;
                if (Mont128::ge128(ys_lo, ys_hi, x_lo, x_hi))
                    Mont128::sub128(diff_lo, diff_hi, ys_lo, ys_hi, x_lo, x_hi);
                else
                    Mont128::sub128(diff_lo, diff_hi, x_lo, x_hi, ys_lo, ys_hi);

                // fromMont, then GCD
                uint64_t dd_lo, dd_hi;
                mont.redc(dd_lo, dd_hi, diff_lo, diff_hi, 0, 0);
                d = gcd128_to64(dd_lo, dd_hi, n_lo, n_hi);
            }
        }

        if (d > 1 && (n_hi > 0 || d < n_lo)) return d;
    }
    return 0;  // Failure
}

} // anonymous namespace

// ============================================================================
// findFactor_PollardRho (Brent's improvement)
// Dispatches to the uint64/128-bit Montgomery/Int version according to bit width
// ============================================================================

Int IntPrime::findFactor_PollardRho(const Int& n) {
    if (n.isEven()) {
        return Int(2);
    }
    if (n <= 1) {
        return n;
    }
    if (isProbablePrime(n)) {
        return n;  // n itself is prime
    }

    // Dispatch to the specialized version according to bit width
    auto bits = n.bitLength();
    if (bits <= 63) {
        uint64_t n64 = n.toUInt64();
        uint64_t f = pollardRho64(n64);
        if (f > 1 && f < n64) return Int(f);
        return n;
    }
    if (bits <= 127) {
        uint64_t n_lo = n.word(0);
        uint64_t n_hi = n.word(1);
        uint64_t f = pollardRho128(n_lo, n_hi);
        if (f > 1) return Int(f);
        return n;
    }

    // 3 limbs or more: the conventional Int-based implementation
    // Brent's rho algorithm
    // Try multiple values of c
    for (int c_val = 1; c_val <= 20; ++c_val) {
        Int c(c_val);
        Int x(2);
        Int y(2);
        Int d(1);

        // f(x) = x² + c mod n
        auto f = [&](const Int& val) -> Int {
            return (val * val + c) % n;
        };

        // Brent's cycle detection with batch GCD
        Int ys, q(1);
        uint64_t m = 128;
        uint64_t r = 1;

        while (d.isOne()) {
            x = y;
            for (uint64_t i = 0; i < r; ++i) {
                y = f(y);
            }

            uint64_t k = 0;
            while (k < r && d.isOne()) {
                ys = y;
                uint64_t batch = (r - k < m) ? (r - k) : m;
                for (uint64_t i = 0; i < batch; ++i) {
                    y = f(y);
                    Int diff = x - y;
                    if (diff.isNegative()) diff = -diff;
                    q = (q * diff) % n;
                }
                d = IntGCD::gcd(q, n);
                k += batch;
            }
            r *= 2;

            // Prevent infinite loops
            if (r > 1000000) break;
        }

        if (d == n) {
            // Backtrack: take the GCD one element at a time
            d = 1;
            while (d.isOne()) {
                ys = f(ys);
                Int diff = x - ys;
                if (diff.isNegative()) diff = -diff;
                d = IntGCD::gcd(diff, n);
            }
        }

        if (d > 1 && d < n) {
            return d;
        }
    }

    return n;  // No factor found
}

// ============================================================================
// findFactor_PollardP1 — Pollard p-1 method
// ============================================================================

Int IntPrime::findFactor_PollardP1(const Int& n, uint64_t B1) {
    if (n.isEven()) return Int(2);
    if (n <= 1) return n;
    if (isProbablePrime(n)) return n;

    // Generate the small-prime list (sieve)
    std::vector<uint64_t> primes;
    {
        std::vector<bool> sieve(B1 + 1, true);
        sieve[0] = sieve[1] = false;
        for (uint64_t i = 2; i * i <= B1; ++i)
            if (sieve[i])
                for (uint64_t j = i * i; j <= B1; j += i)
                    sieve[j] = false;
        for (uint64_t i = 2; i <= B1; ++i)
            if (sieve[i])
                primes.push_back(i);
    }

    // Run Stage 1 with base a=2
    // Incremental GCD checks: check the GCD per prime for early detection
    // (if p-1 of both factors is smooth, the final GCD becomes n)
    Int a(2);
    for (size_t i = 0; i < primes.size(); ++i) {
        uint64_t p = primes[i];
        uint64_t pe = p;
        while (pe <= B1 / p) pe *= p;
        a = IntModular::powerMod(a, Int(pe), n);

        // GCD check at regular intervals (every-time would be heavy, so every 10 primes)
        if ((i + 1) % 10 == 0 || i + 1 == primes.size()) {
            Int a1 = a - 1;
            if (a1.isZero()) break;  // a ≡ 1 (mod n) → no further primes needed
            if (a1.isNegative()) a1 = -a1;
            Int g = IntGCD::gcd(a1, n);
            if (g > 1 && g < n) return g;
            if (g == n) break;  // p-1 of all factors is smooth → retry with a different base
        }
    }

    // If base a=2 failed, retry with other bases
    for (int a_val : {3, 5, 7, 11, 13}) {
        Int a2(a_val);
        for (size_t i = 0; i < primes.size(); ++i) {
            uint64_t p = primes[i];
            uint64_t pe = p;
            while (pe <= B1 / p) pe *= p;
            a2 = IntModular::powerMod(a2, Int(pe), n);

            if ((i + 1) % 10 == 0 || i + 1 == primes.size()) {
                Int a1 = a2 - 1;
                if (a1.isZero()) break;
                if (a1.isNegative()) a1 = -a1;
                Int g = IntGCD::gcd(a1, n);
                if (g > 1 && g < n) return g;
                if (g == n) break;
            }
        }
    }

    return n;  // No factor found
}

// ============================================================================
// findFactor_WilliamsP1 — Williams p+1 method
// Finds a prime factor p where p+1 is B-smooth (the dual of the p-1 method)
// Uses the Lucas sequence V_k(a) mod n
// ============================================================================

namespace {

// Binary-ladder computation of the Lucas sequence V_k(a) mod n
// V_0 = 2, V_1 = a, V_{k+1} = a*V_k - V_{k-1}
// Properties: V_{2k} = V_k² - 2, V_{2k+1} = V_k * V_{k+1} - a
// Keep the pair (V_k, V_{k+1}) and scan the bits of k from the most significant
std::pair<Int, Int> lucasVBinaryLadder(const Int& a, const Int& k, const Int& n) {
    if (k.isZero()) return {Int(2), a};
    if (k == 1) return {a, (a * a - 2) % n};

    Int V0 = a; // V_1
    Int V1 = (a * a - 2) % n; // V_2
    if (V1.isNegative()) V1 += n;

    int bits = static_cast<int>(k.bitLength());
    for (int i = bits - 2; i >= 0; --i) {
        if (k.getBit(i)) {
            // (V_k, V_{k+1}) → (V_{2k+1}, V_{2k+2})
            // V_{2k+1} = V_k * V_{k+1} - a
            // V_{2k+2} = V_{k+1}² - 2
            Int new_V0 = (V0 * V1 - a) % n;
            Int new_V1 = (V1 * V1 - 2) % n;
            V0 = new_V0;
            V1 = new_V1;
        } else {
            // (V_k, V_{k+1}) → (V_{2k}, V_{2k+1})
            // V_{2k} = V_k² - 2
            // V_{2k+1} = V_k * V_{k+1} - a
            Int new_V1 = (V0 * V1 - a) % n;
            Int new_V0 = (V0 * V0 - 2) % n;
            V0 = new_V0;
            V1 = new_V1;
        }
        if (V0.isNegative()) V0 += n;
        if (V1.isNegative()) V1 += n;
    }
    return {V0, V1};
}

} // anonymous namespace

Int IntPrime::findFactor_WilliamsP1(const Int& n, uint64_t B1) {
    if (n.isEven()) return Int(2);
    if (n <= 1) return n;
    if (isProbablePrime(n)) return n;

    // Generate the small-prime list
    std::vector<uint64_t> primes;
    {
        std::vector<bool> sieve(B1 + 1, true);
        sieve[0] = sieve[1] = false;
        for (uint64_t i = 2; i * i <= B1; ++i)
            if (sieve[i])
                for (uint64_t j = i * i; j <= B1; j += i)
                    sieve[j] = false;
        for (uint64_t i = 2; i <= B1; ++i)
            if (sieve[i])
                primes.push_back(i);
    }

    // Try multiple bases a
    for (int a_val : {3, 5, 7, 11, 13, 17, 19, 23}) {
        Int a(a_val);

        // Skip if gcd(a²-4, n) is not 1 (degenerate case)
        Int disc = a * a - 4;
        if (disc.isNegative()) disc = -disc;
        Int g0 = IntGCD::gcd(disc, n);
        if (g0 > 1 && g0 < n) return g0;
        if (g0 == n) continue;

        // V = a (= V_1)
        // For each prime power p^e ≤ B1, compute V = V_{p^e}(V)
        Int V = a;
        for (size_t i = 0; i < primes.size(); ++i) {
            uint64_t p = primes[i];
            uint64_t pe = p;
            while (pe <= B1 / p) pe *= p;

            // V = V_{pe}(V) mod n (composition of the Lucas sequence)
            auto [Vk, Vk1] = lucasVBinaryLadder(V, Int(pe), n);
            V = Vk;

            // Incremental GCD check
            if ((i + 1) % 10 == 0 || i + 1 == primes.size()) {
                Int v2 = V - 2;
                if (v2.isNegative()) v2 += n;
                if (v2.isZero()) break;  // V ≡ 2 (mod n)
                Int g = IntGCD::gcd(v2, n);
                if (g > 1 && g < n) return g;
                if (g == n) break;
            }
        }
    }

    return n;
}

// ============================================================================
// findFactor_SQUFOF — Shanks's square-form factorization method
// An O(n^{1/4}) algorithm based on continued-fraction expansion
// ============================================================================

Int IntPrime::findFactor_SQUFOF(const Int& n) {
    if (n.isEven()) return Int(2);
    if (n <= 1) return n;
    if (isProbablePrime(n)) return n;

    // Perfect-square check
    Int sqrtN = IntSqrt::sqrt(n);
    if (sqrtN * sqrtN == n) return sqrtN;

    // Try small multipliers k (k*n must not be a perfect square)
    static const int multipliers[] = {1, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43};

    for (int k : multipliers) {
        Int D = n * k;

        // Skip if D is a perfect square
        Int sqrtD = IntSqrt::sqrt(D);
        if (sqrtD * sqrtD == D) continue;

        Int P_prev = sqrtD;
        Int Q_prev(1);
        Int Q_curr = D - sqrtD * sqrtD;

        if (Q_curr.isZero()) continue;

        // Forward phase: until Q becomes a perfect square at an even step
        Int P_curr = P_prev;
        int max_iter = 2 * static_cast<int>(std::sqrt(std::sqrt(
            n.bitLength() < 64 ? static_cast<double>(n.toUInt64()) : 1e18)));
        if (max_iter < 1000) max_iter = 1000;
        if (max_iter > 1000000) max_iter = 1000000;

        bool found = false;
        Int q_square_root;

        for (int i = 1; i <= max_iter; ++i) {
            // b_i = floor((sqrtD + P_curr) / Q_curr)
            Int b = (sqrtD + P_curr) / Q_curr;

            // P_{i+1} = b * Q_curr - P_curr
            Int P_next = b * Q_curr - P_curr;

            // Q_{i+1} = Q_{prev} + b * (P_curr - P_next)
            Int Q_next = Q_prev + b * (P_curr - P_next);

            Q_prev = Q_curr;
            Q_curr = Q_next;
            P_curr = P_next;

            // Test Q_{i+1} at even indices (when i is odd, i+1 is even)
            if (i % 2 == 1 && Q_curr > 1) {
                Int sq;
                if (IntSqrt::isSquare(Q_curr, &sq)) {
                    q_square_root = sq;
                    found = true;
                    break;
                }
            }
        }

        if (!found) continue;

        // Reverse phase
        Int q = q_square_root;
        // b_0' = floor((sqrtD - P_curr) / q)
        Int b = (sqrtD - P_curr) / q;

        Int P0 = b * q + P_curr;
        Int Q0 = q;
        Int Q1 = (D - P0 * P0) / Q0;

        if (Q1.isZero()) continue;

        Int Pp = P0, Qp = Q0, Qc = Q1;

        for (int i = 0; i < max_iter; ++i) {
            b = (sqrtD + Pp) / Qc;
            Int P_next = b * Qc - Pp;

            if (P_next == Pp) {
                // Converged: factor = gcd(n, P_next) or gcd(n, Qc)
                Int g = IntGCD::gcd(P_next, n);
                if (g > 1 && g < n) return g;
                g = IntGCD::gcd(Qc, n);
                if (g > 1 && g < n) return g;
                break;
            }

            Int Q_next = Qp + b * (Pp - P_next);
            Qp = Qc;
            Qc = Q_next;
            Pp = P_next;
        }
    }

    return n;
}

// ============================================================================
// findFactor_TrialDivision — trial division method
// ============================================================================

Int IntPrime::findFactor_TrialDivision(const Int& n) {
    if (n <= 1) return n;
    if (n.isEven()) return Int(2);

    // Trial division using the SMALL_PRIMES table (3 to 1987)
    for (int i = 0; i < SMALL_PRIMES_COUNT; ++i) {
        Int p(SMALL_PRIMES[i]);
        if (p * p > n) break;  // Stop once √n is exceeded
        if ((n % p).isZero()) return p;
    }

    // Beyond the table: trial division with odd numbers from 1987 onward (up to √n)
    Int d(1987 + 2);
    while (d * d <= n) {
        if ((n % d).isZero()) return d;
        d += 2;
    }

    return n;  // Prime
}

// ============================================================================
// findFactor_RhoMontgomery — Pollard rho method with Montgomery multiplication
// ============================================================================

namespace {

// -m^{-1} mod 2^64 (Newton iteration)
inline uint64_t rho_mont_neg_inv(uint64_t m0) {
    uint64_t x = 1;
    for (int i = 0; i < 6; ++i)
        x *= 2 - m0 * x;
    return static_cast<uint64_t>(0) - x;
}

// FIOS Montgomery multiplication: r = a * b * R^{-1} mod m
// Fuses multiplication and reduction; no external scratch needed
inline void rho_mont_mul(uint64_t* r,
                         const uint64_t* a, const uint64_t* b,
                         const uint64_t* m, size_t n,
                         uint64_t m_inv) {
    // t[0..n+1]: shift accumulator (stack-allocated)
    uint64_t* t = static_cast<uint64_t*>(_alloca((n + 2) * sizeof(uint64_t)));
    std::memset(t, 0, (n + 2) * sizeof(uint64_t));

#if defined(_MSC_VER) && defined(_M_X64)
    for (size_t i = 0; i < n; ++i) {
        uint64_t bi = b[i];
        uint64_t hi_a, lo_a, hi_m, lo_m;

        lo_a = _umul128(a[0], bi, &hi_a);
        uint64_t S = t[0] + lo_a;
        uint64_t Ca = hi_a + (S < t[0] ? 1ULL : 0);
        unsigned Ca_hi = 0;

        uint64_t q = S * m_inv;

        lo_m = _umul128(m[0], q, &hi_m);
        uint64_t tmp = S + lo_m;
        uint64_t Cm = hi_m + (tmp < S ? 1ULL : 0);
        unsigned Cm_hi = 0;

        for (size_t j = 1; j < n; ++j) {
            lo_a = _umul128(a[j], bi, &hi_a);
            uint64_t s1 = t[j] + lo_a;
            unsigned c1 = (s1 < t[j]) ? 1u : 0u;
            S = s1 + Ca;
            c1 += (S < s1) ? 1u : 0u;
            uint64_t new_Ca = hi_a + c1;
            unsigned new_Ca_hi = (new_Ca < c1) ? 1u : 0u;
            new_Ca += Ca_hi;
            new_Ca_hi += (new_Ca < Ca_hi) ? 1u : 0u;

            lo_m = _umul128(m[j], q, &hi_m);
            uint64_t s2 = S + lo_m;
            unsigned c2 = (s2 < S) ? 1u : 0u;
            uint64_t res = s2 + Cm;
            c2 += (res < s2) ? 1u : 0u;
            uint64_t new_Cm = hi_m + c2;
            unsigned new_Cm_hi = (new_Cm < c2) ? 1u : 0u;
            new_Cm += Cm_hi;
            new_Cm_hi += (new_Cm < Cm_hi) ? 1u : 0u;

            t[j - 1] = res;
            Ca = new_Ca;
            Ca_hi = new_Ca_hi;
            Cm = new_Cm;
            Cm_hi = new_Cm_hi;
        }

        uint64_t f = t[n] + Ca;
        unsigned fc = (f < t[n]) ? 1u : 0u;
        f += Cm;
        fc += (f < Cm) ? 1u : 0u;
        t[n - 1] = f;
        t[n] = static_cast<uint64_t>(fc + Ca_hi + Cm_hi) + t[n + 1];
        t[n + 1] = 0;
    }
#else
    // GCC/Clang: __uint128_t
    for (size_t i = 0; i < n; ++i) {
        uint64_t bi = b[i];
        __uint128_t uv = (__uint128_t)a[0] * bi + t[0];
        uint64_t S = (uint64_t)uv;
        __uint128_t Ca = uv >> 64;

        uint64_t q = S * m_inv;

        uv = (__uint128_t)m[0] * q + S;
        __uint128_t Cm = uv >> 64;

        for (size_t j = 1; j < n; ++j) {
            uv = (__uint128_t)a[j] * bi + t[j] + Ca;
            S = (uint64_t)uv;
            Ca = uv >> 64;

            uv = (__uint128_t)m[j] * q + S + Cm;
            t[j - 1] = (uint64_t)uv;
            Cm = uv >> 64;
        }

        uv = (__uint128_t)t[n] + Ca + Cm;
        t[n - 1] = (uint64_t)uv;
        t[n] = (uint64_t)(uv >> 64) + t[n + 1];
        t[n + 1] = 0;
    }
#endif

    // Conditional subtraction
    if (t[n] != 0 || mpn::cmp(t, n, m, n) >= 0) {
        mpn::sub(r, t, n, m, n);
    } else {
        std::memcpy(r, t, n * sizeof(uint64_t));
    }
}

// mpn level: r += 1, and if r >= m then r -= m
inline void rho_add1_mod(uint64_t* r, const uint64_t* m, size_t n) {
    uint64_t carry = mpn::add_1(r, n, 1);
    if (carry || mpn::cmp(r, n, m, n) >= 0)
        mpn::sub(r, r, n, m, n);
}

} // anonymous namespace

Int IntPrime::findFactor_RhoMontgomery(const Int& n) {
    if (n.isEven()) return Int(2);
    if (n <= 1) return n;
    if (isProbablePrime(n)) return n;

    // mpn-level setup
    const size_t nw = n.size();             // Word count
    const uint64_t* mdata = n.data();

    // m_inv = -m^{-1} mod 2^64
    uint64_t m_inv = rho_mont_neg_inv(mdata[0]);

    // Compute R² mod m
    std::vector<uint64_t> R2(nw, 0);
    if (nw == 1) {
        // nw=1: R = 2^64, computed quickly via _udiv128 / UInt128::divmod_fast
        auto [q1, Rmod1] = UInt128::divmod_fast(1ULL, 0ULL, mdata[0]);
        uint64_t hi, lo;
#if defined(_MSC_VER) && defined(_M_X64)
        lo = _umul128(Rmod1, Rmod1, &hi);
#else
        __uint128_t prod = static_cast<__uint128_t>(Rmod1) * Rmod1;
        hi = static_cast<uint64_t>(prod >> 64);
        lo = static_cast<uint64_t>(prod);
#endif
        auto [q2, R2mod1] = UInt128::divmod_fast(hi, lo, mdata[0]);
        R2[0] = R2mod1;
    } else {
        // nw >= 2: R = B^nw, computed via mpn::divide (requires bn >= 2)
        std::vector<uint64_t> R_buf(nw + 1, 0);
        R_buf[nw] = 1;  // R = 2^(64*nw)
        std::vector<uint64_t> q_buf(2, 0);
        std::vector<uint64_t> Rmod(nw, 0);
        size_t div_sz = mpn::divide_scratch_size(nw + 1, nw);
        size_t mul_sz = mpn::multiply_scratch_size(nw, nw);
        std::vector<uint64_t> work(std::max({div_sz, mul_sz, static_cast<size_t>(1)}), 0);

        mpn::divide(q_buf.data(), Rmod.data(), R_buf.data(), nw + 1, mdata, nw, work.data());

        // (R mod m)² mod m
        size_t rmod_n = nw;
        while (rmod_n > 0 && Rmod[rmod_n - 1] == 0) --rmod_n;
        if (rmod_n == 0) rmod_n = 1;

        std::vector<uint64_t> r2_full(2 * nw, 0);
        mpn::multiply(r2_full.data(), Rmod.data(), rmod_n, Rmod.data(), rmod_n, work.data());
        size_t r2fn = 2 * rmod_n;
        while (r2fn > 0 && r2_full[r2fn - 1] == 0) --r2fn;
        if (r2fn == 0) r2fn = 1;

        if (r2fn <= nw) {
            std::memcpy(R2.data(), r2_full.data(), r2fn * sizeof(uint64_t));
        } else {
            size_t div_sz2 = mpn::divide_scratch_size(r2fn, nw);
            if (div_sz2 > work.size()) work.resize(div_sz2, 0);
            std::vector<uint64_t> q2(r2fn - nw + 1, 0);
            mpn::divide(q2.data(), R2.data(), r2_full.data(), r2fn, mdata, nw, work.data());
        }
    }

    // Work buffers (no heap allocation during the loop)
    std::vector<uint64_t> x(nw, 0), y(nw, 0), temp(nw, 0), diff(nw, 0);

    // Convert the initial value 2 to Montgomery form: Mont(2) = 2 * R mod m = FIOS(2_pad, R²)
    {
        std::vector<uint64_t> two_pad(nw, 0);
        two_pad[0] = 2;
        rho_mont_mul(x.data(), two_pad.data(), R2.data(), mdata, nw, m_inv);
    }
    std::memcpy(y.data(), x.data(), nw * sizeof(uint64_t));

    // Floyd's cycle detection: x = f(x), y = f(f(y))
    // f(x) = x² + 1 (Montgomery form)
    constexpr int max_iter = 1000000;
    for (int iter = 0; iter < max_iter; ++iter) {
        // x = x² + 1 mod n
        rho_mont_mul(temp.data(), x.data(), x.data(), mdata, nw, m_inv);
        std::memcpy(x.data(), temp.data(), nw * sizeof(uint64_t));
        rho_add1_mod(x.data(), mdata, nw);

        // y = f(f(y))
        rho_mont_mul(temp.data(), y.data(), y.data(), mdata, nw, m_inv);
        std::memcpy(y.data(), temp.data(), nw * sizeof(uint64_t));
        rho_add1_mod(y.data(), mdata, nw);
        rho_mont_mul(temp.data(), y.data(), y.data(), mdata, nw, m_inv);
        std::memcpy(y.data(), temp.data(), nw * sizeof(uint64_t));
        rho_add1_mod(y.data(), mdata, nw);

        // diff = |x - y|
        int cmp = mpn::cmp(x.data(), nw, y.data(), nw);
        if (cmp == 0) continue;
        if (cmp > 0)
            mpn::sub(diff.data(), x.data(), nw, y.data(), nw);
        else
            mpn::sub(diff.data(), y.data(), nw, x.data(), nw);

        // Normalize the size of diff and convert to Int (Int only for the GCD call)
        size_t dn = nw;
        while (dn > 0 && diff[dn - 1] == 0) --dn;
        if (dn == 0) continue;

        Int diffInt = Int::fromRawWords(std::span<const uint64_t>(diff.data(), dn), 1);
        Int d = IntGCD::gcd(diffInt, n);

        if (d > 1 && d < n) return d;
        if (d == n) break;  // Failure; a retry is needed
    }

    return n;
}

// ============================================================================
// findFactor_Goldbach — factor search based on the Goldbach conjecture
// Algorithm credit: Dr. Roger G. Doss, PhD
//
// For N = p * q (semiprime):
//   1. Estimate sum ≈ p+q by binary search
//   2. If sum² - 4N is a perfect square, p and q can be solved
// ============================================================================

Int IntPrime::findFactor_Goldbach(const Int& n) {
    if (n.isEven()) return Int(2);
    if (n <= 1) return n;
    if (isProbablePrime(n)) return n;

    // Perfect-square check
    Int sqrtN = IntSqrt::sqrt(n);
    if (sqrtN * sqrtN == n) return sqrtN;

    // Quadratic equation: the solutions of t² - sum*t + N = 0 are p and q
    // Discriminant: solvable if sum² - 4N is a perfect square
    auto tryQuadratic = [&](const Int& sum) -> Int {
        Int disc = sum * sum - 4 * n;
        if (disc.isNegative()) return Int(0);

        Int sq;
        if (!IntSqrt::isSquare(disc, &sq)) return Int(0);

        Int p = (sum - sq) / 2;
        Int q = (sum + sq) / 2;

        if (p > 1 && p < n && p * q == n) return p;
        return Int(0);
    };

    // maxGB: for an estimate of sum ≈ p+q, returns the maximum value of t*s
    // t+s = sum, t decreases from sum/2, s increases from sum/2
    auto maxGB = [](const Int& sum) -> Int {
        Int t = sum >> 1;
        Int s = t;
        while (t > 0) {
            t -= 1;
            s += 1;
            // Return t*s for the first (t, s) pair
            // (sangi also performs a primality check, but it is omitted here)
            return t * s;
        }
        return Int(0);
    };

    // BinarySearch: binary search for sum ≈ p+q
    Int emin(2), emax = n + 1;
    if (emin.isOdd()) emin += 1;
    if (emax.isOdd()) emax += 1;

    Int sum, prev_sum;
    int max_iter = 10000;

    for (int iter = 0; iter < max_iter; ++iter) {
        sum = (emin + emax) >> 1;
        if (sum.isOdd()) sum += 1;

        if (sum == prev_sum) break;  // Converged

        // Try whether the quadratic equation is solvable
        Int result = tryQuadratic(sum);
        if (result > 1 && result < n) return result;

        Int tmp = maxGB(sum);
        if (tmp > n) {
            emax = sum - 1;
        } else if (tmp < n) {
            emin = sum + 1;
        } else {
            // tmp == n: solved directly
            result = tryQuadratic(sum);
            if (result > 1 && result < n) return result;
        }

        prev_sum = sum;
    }

    // After convergence, do a linear search of the neighborhood
    if (sum.isOdd()) sum += 1;
    Int search_end = sum + 10000;
    for (Int s = sum; s < search_end; s += 2) {
        Int result = tryQuadratic(s);
        if (result > 1 && result < n) return result;
    }
    // Reverse direction
    search_end = sum - 10000;
    if (search_end < 2) search_end = Int(2);
    for (Int s = sum; s > search_end; s -= 2) {
        Int result = tryQuadratic(s);
        if (result > 1 && result < n) return result;
    }

    return n;
}

// ============================================================================
// findFactor_ECM — elliptic curve method (Lenstra ECM)
// Montgomery form: By² = x³ + Ax² + x (mod n)
// Scalar multiplication of points via the Montgomery ladder
// ============================================================================

// ECM internals: projective coordinates (X : Z) on a Montgomery curve
// The y coordinate is unnecessary (the Montgomery ladder tracks only the x coordinate)
namespace {

struct MontPoint {
    Int x, z;
};

// Montgomery ladder: point addition (differential addition)
// P3 = P1 + P2 where P0 = P1 - P2 (the difference is known)
// X3 = Z0 * ((X1-Z1)(X2+Z2) + (X1+Z1)(X2-Z2))²
// Z3 = X0 * ((X1-Z1)(X2+Z2) - (X1+Z1)(X2-Z2))²
inline MontPoint montAdd(const MontPoint& P1, const MontPoint& P2,
                          const MontPoint& P0, const Int& n) {
    Int u = (P1.x - P1.z) * (P2.x + P2.z);
    Int v = (P1.x + P1.z) * (P2.x - P2.z);
    Int add = u + v;
    Int sub = u - v;
    MontPoint result;
    result.x = (P0.z * add * add) % n;
    result.z = (P0.x * sub * sub) % n;
    if (result.x.isNegative()) result.x += n;
    if (result.z.isNegative()) result.z += n;
    return result;
}

// Montgomery ladder: point doubling
// X2 = (X1+Z1)² * (X1-Z1)²
// Z2 = ((X1+Z1)² - (X1-Z1)²) * ((X1-Z1)² + A24 * ((X1+Z1)² - (X1-Z1)²))
// where A24 = (A+2)/4
inline MontPoint montDouble(const MontPoint& P, const Int& a24, const Int& n) {
    Int sum = P.x + P.z;
    Int diff = P.x - P.z;
    Int sum2 = sum * sum;
    Int diff2 = diff * diff;
    Int delta = sum2 - diff2;  // = 4*X*Z
    MontPoint result;
    result.x = (sum2 * diff2) % n;
    result.z = (delta * (diff2 + a24 * delta)) % n;
    if (result.x.isNegative()) result.x += n;
    if (result.z.isNegative()) result.z += n;
    return result;
}

// Montgomery ladder: scalar multiple kP
MontPoint montMul(const MontPoint& P, const Int& k, const Int& a24, const Int& n) {
    if (k.isZero()) return {Int(0), Int(0)};
    if (k == 1) return P;

    // Binary method (left-to-right)
    // R0 = P, R1 = 2P
    MontPoint R0 = P;
    MontPoint R1 = montDouble(P, a24, n);

    // Scan the bits of k from the most significant
    int bits = static_cast<int>(k.bitLength());
    for (int i = bits - 2; i >= 0; --i) {
        if (k.getBit(i)) {
            R0 = montAdd(R0, R1, P, n);
            R1 = montDouble(R1, a24, n);
        } else {
            R1 = montAdd(R0, R1, P, n);
            R0 = montDouble(R0, a24, n);
        }
    }
    return R0;
}

} // anonymous namespace

// Run ECM on a single curve — called from the thread pool
// Terminate early once found becomes true
namespace {
Int ecmSingleCurve(const Int& n, int nBits, uint64_t sigma,
                    uint64_t B1, uint64_t B2,
                    const std::vector<uint64_t>& primes,
                    const std::atomic<bool>& found) {
    Int s(sigma);
    Int u = (s * s - 5) % n;
    if (u.isNegative()) u += n;
    Int v = (4 * s) % n;

    MontPoint Q;
    Q.x = (u * u % n * u) % n;
    Q.z = (v * v % n * v) % n;

    Int diff_vu = v - u;
    if (diff_vu.isNegative()) diff_vu += n;
    Int diff3 = (diff_vu * diff_vu % n * diff_vu) % n;
    Int t = (3 * u + v) % n;
    Int a24_num = (diff3 * t) % n;
    Int a24_den = (16 * Q.x % n * v) % n;

    Int g = IntGCD::gcd(a24_den.isNegative() ? -a24_den : a24_den, n);
    if (g > 1 && g < n) return g;
    if (g == n) return Int(0);

    Int inv, dummy_y;
    Int gcd_val = IntGCD::extendedGcd(a24_den, n, inv, dummy_y);
    if (gcd_val != 1) {
        if (gcd_val > 1 && gcd_val < n) return gcd_val;
        return Int(0);
    }
    Int a24 = (a24_num * inv) % n;
    if (a24.isNegative()) a24 += n;

    // Stage 1
    {
        Int zProd(1);
        int gcdInterval = (nBits <= 32) ? 5 : (nBits <= 64) ? 15 : 25;
        int stepsInBatch = 0;
        bool degenerate = false;

        for (uint64_t p : primes) {
            if (p > B1) break;
            if (found.load(std::memory_order_relaxed)) return Int(0);

            uint64_t pe = p;
            while (pe <= B1 / p) pe *= p;
            Q = montMul(Q, Int(pe), a24, n);

            if (Q.z.isZero()) { degenerate = true; break; }

            zProd = (zProd * Q.z) % n;
            ++stepsInBatch;

            if (stepsInBatch >= gcdInterval) {
                g = IntGCD::gcd(zProd.isNegative() ? -zProd : zProd, n);
                if (g > 1 && g < n) return g;
                if (g == n) { degenerate = true; break; }
                zProd = 1;
                stepsInBatch = 0;
            }
        }

        if (degenerate) return Int(0);

        if (stepsInBatch > 0) {
            g = IntGCD::gcd(zProd.isNegative() ? -zProd : zProd, n);
            if (g > 1 && g < n) return g;
            if (g == n) return Int(0);
        }
    }

    g = IntGCD::gcd(Q.z.isNegative() ? -Q.z : Q.z, n);
    if (g > 1 && g < n) return g;
    if (g == n || Q.z.isZero()) return Int(0);

    if (found.load(std::memory_order_relaxed)) return Int(0);

    // Stage 2: Baby-step Giant-step (BSGS)
    // On a Montgomery curve, P and -P have the same x coordinate (X:Z), so
    // X_R·Z_S - Z_R·X_S = 0 detects both kD+r and kD-r
    {
        const uint64_t D = 30;  // 2·3·5 — small, efficient baby-step table

        // Coprime residues of D: gcd(r, 30) = 1, 1 ≤ r < D
        // {1, 7, 11, 13, 17, 19, 23, 29} — 8 of them
        std::vector<uint64_t> residues;
        for (uint64_t r = 1; r < D; r += 2) {
            if (r % 3 == 0 || r % 5 == 0) continue;
            residues.push_back(r);
        }

        // Baby step: build S[r] = [r]Q via differential addition
        std::vector<MontPoint> baby(D);
        baby[1] = Q;
        baby[2] = montDouble(Q, a24, n);
        for (uint64_t r = 3; r < D; r++) {
            baby[r] = montAdd(baby[r-1], baby[1], baby[r-2], n);
        }

        // Obtain [D]Q from the baby steps
        MontPoint QD = baby[D - 1];
        QD = montAdd(QD, baby[1], baby[D - 2], n); // [D]Q = [D-1]Q + Q (diff=[D-2]Q)

        // Start of giant step: k_start·D > B1
        uint64_t k_start = (B1 / D) + 1;
        uint64_t k_end = B2 / D + 1;

        // Compute [k_start·D]Q and [(k_start-1)·D]Q
        MontPoint Rprev = montMul(Q, Int((k_start - 1) * D), a24, n);
        MontPoint Rcur = montMul(Q, Int(k_start * D), a24, n);

        Int prod(1);
        int batch_count = 0;

        for (uint64_t k = k_start; k <= k_end; k++) {
            if (found.load(std::memory_order_relaxed)) return Int(0);

            // For each coprime residue r, check kD ± r
            // Montgomery curve: (X_R·Z_S - Z_R·X_S) = 0 ⟺ [kD]Q = ±[r]Q
            // ⟺ [kD+r]Q = O or [kD-r]Q = O
            for (uint64_t r : residues) {
                // Cover both kD + r and kD - r with a single difference
                uint64_t p_plus = k * D + r;
                uint64_t p_minus = (k * D > r) ? k * D - r : 0;
                if ((p_plus > B2) && (p_minus <= B1 || p_minus == 0)) continue;

                const MontPoint& S = baby[r];
                Int diff_xz = (Rcur.x * S.z - Rcur.z * S.x) % n;
                prod = (prod * diff_xz) % n;
                ++batch_count;
            }

            // Batch GCD check
            if (batch_count >= 200) {
                g = IntGCD::gcd(prod.isNegative() ? -prod : prod, n);
                if (g > 1 && g < n) return g;
                if (g == n) break;
                prod = 1;
                batch_count = 0;
            }

            // Giant step: [(k+1)·D]Q = [kD]Q + [D]Q (diff=[(k-1)·D]Q)
            MontPoint Rnext = montAdd(Rcur, QD, Rprev, n);
            Rprev = Rcur;
            Rcur = Rnext;
        }

        if (batch_count > 0) {
            g = IntGCD::gcd(prod.isNegative() ? -prod : prod, n);
            if (g > 1 && g < n) return g;
        }
    }
    return Int(0);
}
} // anonymous namespace

Int IntPrime::findFactor_ECM(const Int& n, int curves, uint64_t B1, uint64_t B2,
                             const std::atomic<bool>* cancel) {
    if (n.isEven()) return Int(2);
    if (n <= 1) return n;
    if (isProbablePrime(n)) return n;

    int nBits = static_cast<int>(n.bitLength());

    if (B1 == 0) {
        if (nBits <= 20) B1 = 50;
        else if (nBits <= 32) B1 = 150;
        else if (nBits <= 48) B1 = 500;
        else if (nBits <= 64) B1 = 2000;
        else if (nBits <= 96) B1 = 10000;
        else B1 = 50000;
    } else {
        if (nBits <= 20 && B1 > 100) B1 = 100;
        else if (nBits <= 32 && B1 > 500) B1 = 500;
        else if (nBits <= 48 && B1 > 2000) B1 = 2000;
    }

    if (B2 == 0) B2 = B1 * 100;
    if (B2 > B1 * 100) B2 = B1 * 100;

    // Generate the small-prime list
    std::vector<uint64_t> primes;
    {
        std::vector<bool> sieve(B2 + 1, true);
        sieve[0] = sieve[1] = false;
        for (uint64_t i = 2; i * i <= B2; ++i)
            if (sieve[i])
                for (uint64_t j = i * i; j <= B2; j += i)
                    sieve[j] = false;
        for (uint64_t i = 2; i <= B2; ++i)
            if (sieve[i])
                primes.push_back(i);
    }

    // Pre-generate the sigma list (the RNG runs sequentially)
    std::mt19937_64 rng(42);
    uint64_t sigmaRange = (nBits <= 20) ? 50ULL :
                          (nBits <= 48) ? 10000ULL : 1000000ULL;
    std::vector<uint64_t> sigmas(curves);
    for (int i = 0; i < curves; i++)
        sigmas[i] = (rng() % sigmaRange) + 6;

    // Run sequentially for small numbers or a small number of curves (avoids parallelization overhead)
    // For nBits ≤ 48 (~14 digits), the cost per curve is light, so sequential is sufficient
    if (curves < 8 || nBits <= 48) {
        std::atomic<bool> found{false};
        for (int i = 0; i < curves; i++) {
            if (cancel && cancel->load(std::memory_order_relaxed)) return n;
            Int result = ecmSingleCurve(n, nBits, sigmas[i], B1, B2, primes, found);
            if (result > 1 && result < n) return result;
        }
        return n;
    }

    // Parallel execution: distribute the curves over the thread pool
    std::atomic<bool> found{false};
    Int result_factor;
    std::mutex result_mutex;

    std::vector<std::future<void>> futures;
    futures.reserve(curves);

    for (int i = 0; i < curves; i++) {
        futures.push_back(sangi::threadPool().submit([&, sigma = sigmas[i]]() {
            if (found.load(std::memory_order_relaxed)) return;
            if (cancel && cancel->load(std::memory_order_relaxed)) {
                found.store(true, std::memory_order_release);
                return;
            }
            Int f = ecmSingleCurve(n, nBits, sigma, B1, B2, primes, found);
            if (f > 1 && f < n) {
                std::lock_guard<std::mutex> lock(result_mutex);
                if (!found.load(std::memory_order_relaxed)) {
                    result_factor = std::move(f);
                    found.store(true, std::memory_order_release);
                }
            }
        }));
    }

    for (auto& fut : futures)
        fut.get();

    if (found.load(std::memory_order_acquire))
        return result_factor;
    return n;
}

// ============================================================================
// primePi — prime counting via the sieve of Eratosthenes
// ============================================================================

int64_t IntPrime::primePi(int64_t n) {
    if (n < 2) return 0;
    if (n == 2) return 1;

    // Bitmap sieve (odd numbers only)
    int64_t half = n / 2;  // Index k of the odd number 2k+1
    std::vector<bool> sieve(static_cast<size_t>(half + 1), true);
    // sieve[k] = true means (2k+1) is prime

    // Sieve up to sqrt(n)
    for (int64_t i = 1; static_cast<int64_t>((2 * i + 1) * (2 * i + 1)) <= n; ++i) {
        if (sieve[static_cast<size_t>(i)]) {
            int64_t p = 2 * i + 1;  // Prime p
            // Mark from p² in steps of p (multiples of odd numbers only)
            for (int64_t j = (p * p) / 2; j <= half; j += p) {
                sieve[static_cast<size_t>(j)] = false;
            }
        }
    }

    // Count the primes (including 2)
    int64_t count = 1;  // Including 2
    for (int64_t i = 1; 2 * i + 1 <= n; ++i) {
        if (sieve[static_cast<size_t>(i)]) {
            ++count;
        }
    }

    return count;
}

} // namespace sangi