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

// IntPrime.cpp
// 素数判定の実装

#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 calx {

// ============================================================================
// 1リム・バッチ GCD 用テーブルの構築（初回呼び出し時に1度だけ実行）
// ============================================================================

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) {
                // 溢れるので新しいバッチを開始
                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);
        }
        return b;
    }();
    return batches;
}

// ============================================================================
// isDivisibleBySmallPrime — 1リム・バッチ GCD 版
// n が奇数であることが前提（呼び出し元で偶数チェック済み）
// ============================================================================

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

    for (const auto& batch : batches) {
        // N mod batch.product （1リム除算、高速）
        Int r = n % Int(batch.product);
        if (r.isZero()) {
            // batch.product 全体が N を割る → 何かの小素数で割り切れる
            return true;
        }

        // r と batch.product の GCD で因子の有無を判定
        uint64_t rv = r.toUInt64();
        uint64_t g = std::gcd(rv, batch.product);
        if (g > 1) {
            // n 自身が小素数テーブル内の素数かもしれない
            // その場合は "割り切れる" とは言わない
            for (int j = 0; j < batch.count; j++) {
                int p = SMALL_PRIMES[batch.start_index + j];
                if (n == Int(p)) {
                    return false;  // n 自身が素数
                }
            }
            return true;  // n は小素数で割り切れる合成数
        }
    }
    return false;
}

// ============================================================================
// extractSmallPrimeFactors — 1リム・バッチ GCD で小素数因子を一括抽出
// remaining: 入力は奇数（素因数2は事前に除去済み）、出力は小素数因子を除去した値
// factors: 見つかった因子を {素数, 指数} で追加
// ============================================================================

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

    const auto& batches = getBatches();

    // 大きい素数のバッチから先に処理する
    // 大きい素因数で先に割ることで remaining の桁数が早く減り、
    // 後続バッチの N mod batch_product が高速になる
    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 （1リム除算）
        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;  // このバッチの素数はどれも N を割らない

        // g > 1 → このバッチ内のどの素数が remaining を割るか個別チェック
        // 大きい素因数から先に割ることで remaining の桁数を早く減らす
        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;  // この素数は N を割らない

            // p で割り切れるだけ割る
            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 素数判定
// ============================================================================

bool IntPrime::isMillerRabinPrime(const Int& n, int k) {
    // 1. 特殊状態のチェック
    if (n.isNaN() || n.isInfinite()) {
        return false;
    }

    // 2. 2 未満は素数ではない
    if (n < 2) {
        return false;
    }

    // 3. 2 は素数
    if (n == 2) {
        return true;
    }

    // 4. 偶数は素数ではない（2以外）
    if (!n.getBit(0)) {
        return false;
    }

    // 5. 小素数テーブル内の値なら直接判定
    //    SMALL_PRIMES の最大値は 1987
    if (n <= 1987) {
        for (int i = 0; i < SMALL_PRIMES_COUNT; i++) {
            if (n == Int(SMALL_PRIMES[i])) {
                return true;
            }
        }
        // 1987 以下で小素数リストにない奇数 → 合成数
        return false;
    }

    // 6. バッチ GCD で小素数による割り切れチェック
    if (isDivisibleBySmallPrime(n)) {
        return false;
    }

    // 7. n-1 = 2^s × d と分解（d は奇数）
    Int n_minus_1 = n - 1;
    Int d = n_minus_1;
    int s = static_cast<int>(d.countTrailingZeros());
    d = d >> s;

    // 8. k 回の Miller-Rabin テスト
    std::random_device rd;
    std::mt19937_64 gen(rd());

    // 底として小素数を使い、足りなければランダム
    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);
        }
        // 退化した witness (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 素数判定
// ============================================================================

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) {
    // n が偶数なら 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;  // 完全平方数
    }
    a += 1;

    // b² = a² - n を計算し、b² が完全平方数になるまで a をインクリメント
    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;  // 因数が見つからなかった
}

// ============================================================================
// Pollard rho — uint64 専用高速版 (n ≤ 63 bits)
// Int 構築を完全排除し、128-bit 乗算 + binary GCD で高速化
// ============================================================================

namespace {

// 64-bit 剰余乗算: (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 専用、バッチ 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 なので溢れない (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) {
            // バックトラック
            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 専用版 (64 < n ≤ 127 bits)
// ============================================================================

struct Mont128 {
    uint64_t n_lo, n_hi;     // n (128-bit)
    uint64_t inv_lo;         // -n^{-1} mod 2^64 (Montgomery 逆元の下位)

    // 128-bit 加算: a + b, 結果が n 以上なら 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 減算: a - b (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 比較: 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 方式: T に m_i * n を加算して下位リムを消去
    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  (最大 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 = (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;  // 新しい t1 (次の 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;
        }

        // 条件付き減算: result >= n なら 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 乗算: 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 結果 (GCD が 64-bit に収まる前提)
uint64_t gcd128_to64(uint64_t a_lo, uint64_t a_hi,
                     uint64_t b_lo, uint64_t b_hi) {
    // 簡易版: Int の GCD を使う
    // rho のバッチ GCD では product を n で剰余した 128-bit 値と n の GCD
    // 結果は因数なので 64-bit 以下
    if (a_hi == 0 && a_lo == 0) {
        if (b_hi == 0) return b_lo;
        // b が大きい → Int で
    }
    if (b_hi == 0 && b_lo == 0) {
        if (a_hi == 0) return a_lo;
    }
    // 両方 64-bit なら高速パス
    if (a_hi == 0 && b_hi == 0) return gcd64(a_lo, b_lo);

    // 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 版
// n は 2 リム (64 < bitLength <= 127, 奇数)
uint64_t pollardRho128(uint64_t n_lo, uint64_t n_hi) {
    // Montgomery セットアップ
    Mont128 mont;
    mont.n_lo = n_lo;
    mont.n_hi = n_hi;

    // -n^{-1} mod 2^64 (Newton 法: x_{i+1} = x_i * (2 - n * x_i))
    // 初期値 x_0 = 1 (奇数 n に対して n*1 ≡ 1 mod 2)
    // 6 反復で 64-bit 精度に収束
    uint64_t inv = 1;
    for (int i = 0; i < 6; ++i)
        inv = inv * (2 - n_lo * inv);
    // inv = n^{-1} mod 2^64, 否定して -n^{-1} mod 2^64
    mont.inv_lo = ~inv + 1;  // = -inv mod 2^64

    // R = 2^128 mod n, R² mod n を計算 (通常の除算で初期化)
    // 簡易: Int 経由で R² mod n を計算
    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 乗算)
    // fromMont(x) = x * 1 (REDC)

    for (uint64_t c_val = 1; c_val <= 20; ++c_val) {
        // c, x, y を Montgomery 形式に変換
        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| (Montgomery 形式のまま差を取る)
                    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 乗算)
                    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) — q を通常形式に戻して 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) {
            // d == n の場合: バックトラック (n_hi != 0 のときは単純比較不可)
            goto backtrack;
        }
        if (n_hi != 0) {
            // d は 64-bit だが n は 128-bit → d == n にはならない
            // ただし d が n_lo と同じで n_hi が 0 でない場合は d < n 確定
        }
        // 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 して 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;  // 失敗
}

} // anonymous namespace

// ============================================================================
// findFactor_PollardRho (Brent's improvement)
// ビット幅に応じて uint64/128-bit Montgomery/Int 版をディスパッチ
// ============================================================================

Int IntPrime::findFactor_PollardRho(const Int& n) {
    if (n.isEven()) {
        return Int(2);
    }
    if (n <= 1) {
        return n;
    }
    if (isProbablePrime(n)) {
        return n;  // n 自体が素数
    }

    // ビット幅に応じて専用版にディスパッチ
    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 リム以上: 従来の Int ベース実装
    // Brent's rho algorithm
    // 複数の 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;

            // 無限ループ防止
            if (r > 1000000) break;
        }

        if (d == n) {
            // バックトラック: 1要素ずつ GCD を取る
            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;  // 因数が見つからなかった
}

// ============================================================================
// findFactor_PollardP1 — Pollard p-1 法
// 移植元: lib++20 Factor_P1 (改善版)
// ============================================================================

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;

    // 小素数リスト生成 (篩)
    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);
    }

    // 底 a=2 で Stage 1 実行
    // 段階的 GCD チェック: 素数ごとに GCD を確認して早期検出
    // (両方の因数の p-1 が smooth だと最後の GCD が 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 チェック (毎回だと重いので 10 素数ごと)
        if ((i + 1) % 10 == 0 || i + 1 == primes.size()) {
            Int a1 = a - 1;
            if (a1.isZero()) break;  // a ≡ 1 (mod n) → これ以上の素数は不要
            if (a1.isNegative()) a1 = -a1;
            Int g = IntGCD::gcd(a1, n);
            if (g > 1 && g < n) return g;
            if (g == n) break;  // 全因数の p-1 が smooth → 別の底で再試行
        }
    }

    // 底 a=2 で失敗した場合、他の底で再試行
    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;  // 因数が見つからなかった
}

// ============================================================================
// findFactor_WilliamsP1 — Williams p+1 法
// p+1 が B-smooth な素因数 p を見つける (p-1 法の双対)
// Lucas 列 V_k(a) mod n を使用
// ============================================================================

namespace {

// Lucas 列 V_k(a) mod n の binary ladder 計算
// V_0 = 2, V_1 = a, V_{k+1} = a*V_k - V_{k-1}
// 性質: V_{2k} = V_k² - 2, V_{2k+1} = V_k * V_{k+1} - a
// (V_k, V_{k+1}) のペアを保持し、k のビットを上位から走査
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;

    // 小素数リスト生成
    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);
    }

    // 複数の底 a で試行
    for (int a_val : {3, 5, 7, 11, 13, 17, 19, 23}) {
        Int a(a_val);

        // gcd(a²-4, n) が 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)
        // 各素数冪 p^e ≤ B1 について 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 (Lucas 列の合成)
            auto [Vk, Vk1] = lucasVBinaryLadder(V, Int(pe), n);
            V = Vk;

            // 段階的 GCD チェック
            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 の平方形式因数分解法
// 連分数展開に基づく O(n^{1/4}) アルゴリズム
// ============================================================================

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

    // 完全平方数チェック
    Int sqrtN = IntSqrt::sqrt(n);
    if (sqrtN * sqrtN == n) return sqrtN;

    // 小さい乗数 k で試行 (k*n が完全平方数にならないこと)
    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;

        // D が完全平方数ならスキップ
        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;

        // 前進フェーズ: 偶数ステップで Q が完全平方数になるまで
        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;

            // Q_{i+1} を偶数インデックスで判定 (i が奇数のとき i+1 は偶数)
            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;

        // 逆行フェーズ
        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) {
                // 収束: factor = gcd(n, P_next) または 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 — 試し割り法
// 移植元: lib++20 Factor_Try
// ============================================================================

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

    // SMALL_PRIMES テーブル (3〜1987) で試し割り
    for (int i = 0; i < SMALL_PRIMES_COUNT; ++i) {
        Int p(SMALL_PRIMES[i]);
        if (p * p > n) break;  // √n を超えたら打ち切り
        if ((n % p).isZero()) return p;
    }

    // テーブル外: 1987 以降の奇数で試し割り (√n まで)
    Int d(1987 + 2);
    while (d * d <= n) {
        if ((n % d).isZero()) return d;
        d += 2;
    }

    return n;  // 素数
}

// ============================================================================
// findFactor_RhoMontgomery — Montgomery 乗算版 Pollard rho 法
// 移植元: lib++20 Factor_Rho_Montgomery
// ============================================================================

namespace {

// Montgomery 剰余: T * R^{-1} mod N
// T < R*N を仮定
// REDC(T) = (T + (T * Ninv mod R) * N) / R
// ただし結果が N 以上なら N を引く
Int montRedc(const Int& T, const Int& N, const Int& Ninv, int pos, const Int& R_mask) {
    // m = (T mod R) * Ninv mod R
    // R = 2^pos なので mod R は下位 pos ビットのマスク (R_mask = R - 1 を事前計算)
    Int t_low = T & R_mask;
    Int m = (t_low * Ninv) & R_mask;

    // result = (T + m * N) >> pos
    Int result = (T + m * N) >> pos;

    if (result >= N) result -= N;
    if (result.isNegative()) result += N;
    return result;
}

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

    // R = 2^pos > n (n のワード数の整数倍ビット)
    int nbits = static_cast<int>(n.bitLength());
    int pos = ((nbits + 63) / 64) * 64;  // 64 ビットの倍数に切り上げ
    Int R = Int(1) << pos;

    // R² mod N (Montgomery 表現への変換に使用)
    Int R2 = (R * R) % n;

    // 拡張ユークリッドで Ninv を求める: R * Rinv - N * Ninv = 1
    // つまり N * Ninv ≡ -1 (mod R)
    Int Rinv, Ninv;
    IntGCD::extendedGcd(R, n, Rinv, Ninv);
    // R*Rinv + N*(-Ninv) = 1 → N*Ninv ≡ -1 (mod R)
    // extendedGcd(R, -N, Rinv, Ninv) としたいが、
    // R*Rinv + n*Ninv = gcd = 1 なので
    // Ninv = -Ninv (n の係数) が N*Ninv ≡ -1 mod R の解
    Ninv = -Ninv;
    Int R_mask = R - 1;
    Ninv = Ninv & R_mask;
    if (Ninv.isNegative()) Ninv += R;

    // Montgomery 表現への変換: aR mod N = REDC(a * R²)
    // x = 2 in Montgomery form
    Int x = montRedc(2 * R2, n, Ninv, pos, R_mask);
    Int y = x;
    Int d;

    // Floyd のサイクル検出: x = f(x), y = f(f(y))
    // f(x) = x² + 1 (Montgomery 乗算)
    int max_iter = 1000000;
    for (int iter = 0; iter < max_iter; ++iter) {
        // x = x² + 1 mod n (Montgomery)
        x = montRedc(x * x, n, Ninv, pos, R_mask);
        x += 1;
        if (x >= n) x -= n;

        // y = f(f(y))
        y = montRedc(y * y, n, Ninv, pos, R_mask);
        y += 1;
        if (y >= n) y -= n;
        y = montRedc(y * y, n, Ninv, pos, R_mask);
        y += 1;
        if (y >= n) y -= n;

        // d = gcd(|x - y|, n)
        Int diff = x - y;
        if (diff.isNegative()) diff = -diff;
        if (diff.isZero()) continue;

        d = IntGCD::gcd(diff, n);

        if (d > 1 && d < n) return d;
        if (d == n) break;  // 失敗、再試行が必要
    }

    return n;
}

// ============================================================================
// findFactor_Goldbach — Goldbach 予想に基づく因数探索
// 移植元: lib++20 Factor_Goldbach (Dr. Roger G. Doss, PhD)
//
// N = p * q (半素数) に対し:
//   1. sum ≈ p+q を二分探索で推定
//   2. sum² - 4N が完全平方数なら p,q を解ける
// ============================================================================

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

    // 完全平方数チェック
    Int sqrtN = IntSqrt::sqrt(n);
    if (sqrtN * sqrtN == n) return sqrtN;

    // 二次方程式: t² - sum*t + N = 0 の解が p と q
    // 判別式: sum² - 4N が完全平方数なら解ける
    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: sum ≈ p+q の推定値に対し、t*s の最大値を返す
    // t+s = sum, t は sum/2 から減少、s は sum/2 から増加
    auto maxGB = [](const Int& sum) -> Int {
        Int t = sum >> 1;
        Int s = t;
        while (t > 0) {
            t -= 1;
            s += 1;
            // 最初の (t, s) ペアで t*s を返す
            // (lib++20 では素数チェックも行うが、ここでは省略)
            return t * s;
        }
        return Int(0);
    };

    // BinarySearch: 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;  // 収束

        // 二次方程式で解けるか試す
        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: 直接解けた
            result = tryQuadratic(sum);
            if (result > 1 && result < n) return result;
        }

        prev_sum = sum;
    }

    // 収束後、周辺を線形探索
    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;
    }
    // 逆方向
    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 — 楕円曲線法 (Lenstra ECM)
// Montgomery 形式: By² = x³ + Ax² + x (mod n)
// Montgomery ladder による点のスカラー倍算
// ============================================================================

// ECM 内部: Montgomery 曲線上の射影座標 (X : Z)
// y 座標は不要（Montgomery ladder では x 座標のみ追跡）
namespace {

struct MontPoint {
    Int x, z;
};

// Montgomery ladder: 点の加算 (differential addition)
// P3 = P1 + P2 where P0 = P1 - P2 (差分が既知)
// 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: 点の2倍
// 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: スカラー倍 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);

    // k のビット列を上位から走査
    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

// ECM 単一曲線の実行 — スレッドプールから呼び出される
// found が 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)
    // Montgomery 曲線では P と -P が同じ x 座標 (X:Z) を持つため
    // X_R·Z_S - Z_R·X_S = 0 で kD+r と kD-r の両方を検出
    {
        const uint64_t D = 30;  // 2·3·5 — baby step テーブルが小さく効率的

        // D の coprime 残基: gcd(r, 30) = 1, 1 ≤ r < D
        // {1, 7, 11, 13, 17, 19, 23, 29} — 8 個
        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: S[r] = [r]Q を差分加算で構築
        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);
        }

        // [D]Q を baby step から取得
        MontPoint QD = baby[D - 1];
        QD = montAdd(QD, baby[1], baby[D - 2], n); // [D]Q = [D-1]Q + Q (diff=[D-2]Q)

        // Giant step 開始: k_start·D > B1
        uint64_t k_start = (B1 / D) + 1;
        uint64_t k_end = B2 / D + 1;

        // [k_start·D]Q と [(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);

            // 各 coprime 残基 r について kD ± r をチェック
            // Montgomery curve: (X_R·Z_S - Z_R·X_S) = 0 ⟺ [kD]Q = ±[r]Q
            // ⟺ [kD+r]Q = O または [kD-r]Q = O
            for (uint64_t r : residues) {
                // kD + r と kD - r の両方を単一の差分でカバー
                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;
            }

            // バッチ GCD チェック
            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;

    // 小素数リスト生成
    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);
    }

    // sigma リストを事前生成 (RNG は逐次実行)
    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;

    // 小さい数 or 曲線数が少ない場合は逐次実行（並列化オーバーヘッド回避）
    // nBits ≤ 48 (~14桁) では 1 曲線あたりのコストが軽く、逐次で十分
    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;
    }

    // 並列実行: 曲線をスレッドプールで分散
    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(calx::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 — エラトステネスの篩による素数計数
// ============================================================================

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

    // ビットマップ篩 (奇数のみ)
    int64_t half = n / 2;  // 奇数 2k+1 のインデックス k
    std::vector<bool> sieve(static_cast<size_t>(half + 1), true);
    // sieve[k] = true means (2k+1) is prime

    // 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;  // 素数 p
            // p² から p ずつ (奇数の倍数のみ) マーク
            for (int64_t j = (p * p) / 2; j <= half; j += p) {
                sieve[static_cast<size_t>(j)] = false;
            }
        }
    }

    // 素数を数える (2 を含む)
    int64_t count = 1;  // 2 を含む
    for (int64_t i = 1; 2 * i + 1 <= n; ++i) {
        if (sieve[static_cast<size_t>(i)]) {
            ++count;
        }
    }

    return count;
}

} // namespace calx