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

// IntPrimeGNFS.cpp
// GNFS (General Number Field Sieve) 簡易実装
// 教育目的の基本版: base-m polynomial, line sieve, Hensel algebraic sqrt

#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/Int/IntNumberTheory.hpp"
#include "math/core/mp/Int/BlockLanczos.hpp"

#include <vector>
#include <cmath>
#include <cstdint>
#include <algorithm>
#include <numeric>
#include <random>
#include <unordered_map>
#include <unordered_set>
#include <queue>
#include <fstream>
#include <sstream>
#include <thread>
#include <mutex>

namespace calx {
namespace {

// ============================================================
// uint64 modular helpers (p < 2^32 保証)
// ============================================================

uint64_t gnfs_powmod(uint64_t b, uint64_t e, uint64_t m) {
    uint64_t r = 1; b %= m;
    while (e) { if (e & 1) r = r * b % m; b = b * b % m; e >>= 1; }
    return r;
}

uint64_t gnfs_modinv(uint64_t a, uint64_t p) {
    int64_t r0 = (int64_t)(a % p), r1 = (int64_t)p, s0 = 1, s1 = 0;
    while (r1) { int64_t q = r0/r1, t; t=r1; r1=r0-q*r1; r0=t; t=s1; s1=s0-q*s1; s0=t; }
    return s0 < 0 ? (uint64_t)(s0 + p) : (uint64_t)s0;
}

// ============================================================
// 64-bit modular arithmetic (for cofactor analysis)
// ============================================================

uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t hi;
    uint64_t lo = _umul128(a, b, &hi);
    uint64_t r;
    _udiv128(hi, lo, m, &r);
    return r;
}

// Deterministic Miller-Rabin for 64-bit
// Witnesses {2,3,5,7,11,13,17,19,23,29,31,37} are sufficient for n < 2^64
bool isPrime64(uint64_t n) {
    if (n < 2) return false;
    if (n < 4) return true;
    if (n % 2 == 0) return false;
    uint64_t d = n - 1;
    int r = 0;
    while (d % 2 == 0) { d /= 2; r++; }
    for (uint64_t a : {2,3,5,7,11,13,17,19,23,29,31,37}) {
        if (a >= n) continue;
        uint64_t x = 1, base = a % n, e = d;
        while (e) { if (e & 1) x = mulmod64(x, base, n); base = mulmod64(base, base, n); e >>= 1; }
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int i = 0; i < r - 1; i++) {
            x = mulmod64(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

// Pollard's rho for 64-bit composite
uint64_t pollardRho64(uint64_t n) {
    if (n % 2 == 0) return 2;
    for (uint64_t c = 1; c < 100; c++) {
        uint64_t x = 2, y = 2, d = 1;
        while (d == 1) {
            x = (mulmod64(x, x, n) + c) % n;
            y = (mulmod64(y, y, n) + c) % n;
            y = (mulmod64(y, y, n) + c) % n;
            d = std::gcd(x > y ? x - y : y - x, n);
        }
        if (d != n) return d;
    }
    return 0;
}

uint64_t gnfs_sqrtmod(uint64_t n, uint64_t p) {
    n %= p; if (n == 0) return 0; if (p == 2) return n & 1;
    if (p % 4 == 3) return gnfs_powmod(n, (p+1)/4, p);
    uint64_t Q = p - 1; int S = 0;
    while (Q % 2 == 0) { Q /= 2; S++; }
    uint64_t z = 2;
    while (gnfs_powmod(z, (p-1)/2, p) != p-1) z++;
    uint64_t M = S, c = gnfs_powmod(z, Q, p);
    uint64_t t = gnfs_powmod(n, Q, p), R = gnfs_powmod(n, (Q+1)/2, p);
    while (true) {
        if (t == 1) return R;
        uint64_t tmp = t; int i = 0;
        while (tmp != 1) { tmp = tmp * tmp % p; i++; }
        uint64_t b = c;
        for (int j = 0; j < (int)M - i - 1; j++) b = b * b % p;
        M = i; c = b * b % p; t = t * c % p; R = R * b % p;
    }
}

// ============================================================
// Polynomial in Z[x]/(f(x)) mod M
// f は monic (f[d]=1), poly は d 個の係数 (degree < d)
// ============================================================

using Poly = std::vector<Int>;

Poly polyMul(const Poly& a, const Poly& b, const Poly& f, int d, const Int& M) {
    std::vector<Int> prod(2*d - 1, Int(0));
    for (int i = 0; i < d; i++)
        for (int j = 0; j < d; j++)
            prod[i+j] = (prod[i+j] + a[i] * b[j]) % M;
    for (int i = 2*d-2; i >= d; i--) {
        if (prod[i].isZero()) continue;
        for (int j = 0; j < d; j++)
            prod[i-d+j] = (prod[i-d+j] - prod[i] * f[j]) % M;
    }
    Poly r(d);
    for (int i = 0; i < d; i++) { r[i] = prod[i] % M; if (r[i].isNegative()) r[i] += M; }
    return r;
}

// GF(2) Gaussian elimination は gf2_linalg::findNullSpaceGF2 (BlockLanczos.hpp) に統合

// ============================================================
// GNFS パラメータ
// ============================================================

struct GNFSParams {
    int degree, rfb_bound, sieve_a, sieve_b, thresh_adj;
};

GNFSParams gnfsParams(int digits) {
    if (digits <= 30) return {3,  3000,  20000,  100,  18};
    if (digits <= 40) return {3,  8000,  50000,  200,  20};
    if (digits <= 50) return {3, 15000,  80000,  400,  22};
    if (digits <= 60) return {4, 30000, 120000,  600,  24};
    return                   {5, 60000, 200000, 1000,  26};
}

// ============================================================
// Base-m polynomial selection
// n = f(m), f monic degree d, m ≈ n^{1/d}
// ============================================================

// Forward declarations
static std::vector<uint64_t> findRootsModP(const Poly& f_int, int d, uint64_t p);

struct GNFSPoly { Poly f; Int m; int d; double score; double skewness = 1.0; };

// ============================================================
// 多項式スコアリング (skewness + root properties)
// ============================================================

// Alpha score: 小素数での根の個数を反映 (Murphy's alpha)
// alpha が小さい (負方向に大きい) ほど smooth 確率が高い
// 拡張版: 200 以下の素数を使用し、正確な寄与を計算
double polyAlpha(const Poly& f, int d) {
    double alpha = 0.0;
    // 200 以下の素数 (46 個) でスコアリング
    uint64_t primes[] = {
        2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
        73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,
        151,157,163,167,173,179,181,191,193,197,199
    };
    for (uint64_t p : primes) {
        auto roots = findRootsModP(f, d, p);
        int rp = static_cast<int>(roots.size());
        // Murphy's alpha: Σ (r_p/(p-1) - 1/p) · log(p)
        // ≈ smooth probability contribution from prime p
        double contrib = static_cast<double>(rp) / static_cast<double>(p - 1)
                       - 1.0 / static_cast<double>(p);
        alpha -= contrib * std::log(static_cast<double>(p));
    }
    return alpha;
}

// Skewness: 最適な s を L2 ノルム最小化で決定
// 初期推定 s₀ = (|c₀|/|c_d|)^{1/d}, 黄金分割探索で精密化
double polySkewness(const Poly& f, int d) {
    double c0 = std::abs(f[0].toDouble());
    double cd = std::abs(f[d].toDouble());
    if (cd < 1.0 || c0 < 1.0) return 1.0;
    double s0 = std::pow(c0 / cd, 1.0 / d);

    // L2 norm as function of skewness
    auto l2norm = [&](double s) {
        double norm2 = 0.0;
        double sp = 1.0;
        for (int i = 0; i <= d; i++) {
            double ci = f[i].toDouble() / sp;
            norm2 += ci * ci;
            sp *= s;
        }
        return norm2;
    };

    // 黄金分割探索で最適 skewness を求める
    constexpr double phi = 1.618033988749895;
    double lo = s0 / 4.0, hi = s0 * 4.0;
    if (lo < 1.0) lo = 1.0;

    for (int iter = 0; iter < 50; iter++) {
        double m1 = hi - (hi - lo) / phi;
        double m2 = lo + (hi - lo) / phi;
        if (l2norm(m1) < l2norm(m2))
            hi = m2;
        else
            lo = m1;
    }
    return (lo + hi) / 2.0;
}

// Murphy E-score: 篩効率の推定 (小さいほど良い)
// E ≈ log₂(skewed_L2_norm) - alpha - root_score
// root_score: 有理側 (m の root properties) も考慮
double polyScore(const Poly& f, int d) {
    double skew = polySkewness(f, d);

    // Skewed L2 norm
    double norm2 = 0.0;
    double sp = 1.0;
    for (int i = 0; i <= d; i++) {
        double ci = f[i].toDouble() / sp;
        norm2 += ci * ci;
        sp *= skew;
    }
    double lognorm = 0.5 * std::log2(std::max(norm2, 1.0));

    // Alpha (root properties of algebraic polynomial)
    double alpha = polyAlpha(f, d);

    // 有理側の寄与: c_d が小さいほど有利 (有理ノルム = |a + bm| が小さい)
    double rat_bonus = 0.0;
    double ad = std::abs(f[d].toDouble());
    if (ad > 0.0 && ad < 100.0) {
        // 先頭係数が小さい → 有理側のスムース確率が上がる
        rat_bonus = std::log2(ad + 1.0) * 0.3;
    }

    return lognorm - alpha + rat_bonus;
}

// n を基数 m で表現し monic 多項式を生成
static bool baseMExpand(const Int& n, const Int& m, int d, Poly& f) {
    f.resize(d + 1);
    Int rem = n;
    for (int i = 0; i < d; i++) { f[i] = rem % m; rem = rem / m; }
    f[d] = rem;
    return (f[d] == 1);
}

GNFSPoly baseMSelect(const Int& n, int d) {
    auto ipow = [](const Int& base, int exp) -> Int {
        Int r(1); for (int i = 0; i < exp; i++) r *= base; return r;
    };

    // Newton's method for d-th root
    double nd = std::log2(std::max(1.0, n.toDouble())) / d;
    Int m0 = Int(static_cast<uint64_t>(std::pow(2.0, nd)));
    if (m0 < 2) m0 = 2;

    // m^d ≤ n < (m+1)^d
    for (int iter = 0; iter < 100; iter++) {
        Int md = ipow(m0, d);
        if (md > n) { m0 -= 1; continue; }
        Int md1 = ipow(m0 + 1, d);
        if (md1 <= n) { m0 += 1; continue; }
        break;
    }

    // Phase 1: m0 ± 200 の範囲でスコア最良の monic 多項式を探索
    GNFSPoly best;
    best.d = d;
    best.score = 1e100;
    bool found = false;

    for (int delta = -200; delta <= 200; delta++) {
        Int m2 = m0 + delta;
        if (m2 < 2) continue;
        Poly f2;
        if (!baseMExpand(n, m2, d, f2)) continue;
        double s = polyScore(f2, d);
        if (s < best.score) {
            best.f = f2;
            best.m = m2;
            best.score = s;
            found = true;
        }
    }

    if (!found) {
        Poly f;
        baseMExpand(n, m0, d, f);
        return {f, m0, d, 0.0};
    }

    // Phase 2: 最良 monic を非 monic (小さい先頭係数) に変形して比較
    // c_d ∈ {2, 3, 5} で n - c_d·m^d の残余を分配
    for (int cd_val = 2; cd_val <= 5; cd_val++) {
        // m' = floor((n / c_d)^{1/d})
        Int n_div = n / Int(cd_val);
        double nd2 = std::log2(std::max(1.0, n_div.toDouble())) / d;
        Int m2_base = Int(static_cast<uint64_t>(std::pow(2.0, nd2)));
        if (m2_base < 2) continue;

        // m'^d · c_d ≤ n の保証
        for (int iter = 0; iter < 100; iter++) {
            Int md = ipow(m2_base, d) * Int(cd_val);
            if (md > n) { m2_base -= 1; continue; }
            Int md1 = ipow(m2_base + 1, d) * Int(cd_val);
            if (md1 <= n) { m2_base += 1; continue; }
            break;
        }

        // m2_base ± 20 で探索
        for (int delta = -20; delta <= 20; delta++) {
            Int m2 = m2_base + delta;
            if (m2 < 2) continue;

            // f(x) = c_d·x^d + c_{d-1}·x^{d-1} + ... + c_0 where f(m2) = n
            Int rem = n;
            Poly f2(d + 1);
            f2[d] = Int(cd_val);
            rem -= Int(cd_val) * ipow(m2, d);
            if (rem.isNegative()) continue;

            bool valid = true;
            for (int i = d - 1; i >= 0; i--) {
                Int mi = ipow(m2, i);
                if (mi.isZero()) { valid = false; break; }
                f2[i] = rem / mi;
                rem -= f2[i] * mi;
                // 負の係数も許容 (|c_i| < m が望ましい)
            }
            if (!valid || !rem.isZero()) continue;

            double s = polyScore(f2, d);
            if (s < best.score) {
                best.f = f2;
                best.m = m2;
                best.score = s;
            }
        }
    }

    return best;
}

// ============================================================
// Kleinjung-style polynomial selection (GNFS-ADV-1)
// ============================================================

// Smooth c_d 列挙: 小素数 {2,3,5,7,11,13} の積で構成される数を生成
// 射影根が多い c_d ほど alpha が改善される
static std::vector<uint64_t> enumerateSmoothCd(int d, uint64_t max_cd) {
    // Smooth 数を BFS 的に生成
    const uint64_t small_primes[] = {2, 3, 5, 7, 11, 13};
    std::vector<uint64_t> smooth;
    smooth.push_back(1);

    // ヒープ的に重複なく生成
    std::vector<uint64_t> heap = {1};
    std::unordered_set<uint64_t> seen;
    seen.insert(1);

    while (!heap.empty()) {
        std::sort(heap.begin(), heap.end());
        uint64_t cur = heap.front();
        heap.erase(heap.begin());

        for (uint64_t p : small_primes) {
            if (cur > max_cd / p) continue;  // overflow 回避
            uint64_t next = cur * p;
            if (next > max_cd) continue;
            if (seen.count(next)) continue;
            seen.insert(next);
            smooth.push_back(next);
            heap.push_back(next);
        }
    }

    std::sort(smooth.begin(), smooth.end());
    return smooth;
}

// 射影根スコア: c_d が各素数 p で何個の射影根を持つか
// c_d の素因数 p は自動的に射影根を 1 つ寄与する
static double projAlphaScore(uint64_t cd) {
    double score = 0.0;
    const uint64_t primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31};
    for (uint64_t p : primes) {
        int div_count = 0;
        uint64_t tmp = cd;
        while (tmp % p == 0) { tmp /= p; div_count++; }
        if (div_count > 0) {
            // 射影根 1 つ分の alpha 改善
            score += std::log(static_cast<double>(p)) / static_cast<double>(p - 1);
        }
    }
    return score;
}

// 2D サイズ縮約 (Lagrange): 格子 {(e1x,e1y), (e2x,e2y)} の短ベクトルを求める
// Int 版 (c_{d-1} が巨大になりうる)
static void latticeReduce2D(Int& e1x, Int& e1y, Int& e2x, Int& e2y) {
    auto norm2d = [](const Int& x, const Int& y) -> Int {
        return x * x + y * y;
    };

    for (int iter = 0; iter < 200; iter++) {
        Int n1 = norm2d(e1x, e1y);
        Int n2 = norm2d(e2x, e2y);
        if (n1 > n2) {
            std::swap(e1x, e2x);
            std::swap(e1y, e2y);
            std::swap(n1, n2);
        }
        if (n1.isZero()) break;
        // dot = e1·e2
        Int dot = e1x * e2x + e1y * e2y;
        // mu = round(dot / n1)
        // 四捨五入: (2·dot + n1) / (2·n1) (dot ≥ 0), (2·dot - n1) / (2·n1) (dot < 0)
        Int mu;
        if (dot.isNegative()) {
            mu = (2 * dot - n1) / (2 * n1);
        } else {
            mu = (2 * dot + n1) / (2 * n1);
        }
        if (mu.isZero()) break;
        e2x -= mu * e1x;
        e2y -= mu * e1y;
    }

    if (norm2d(e1x, e1y) > norm2d(e2x, e2y)) {
        std::swap(e1x, e2x);
        std::swap(e1y, e2y);
    }
}

// 線形回転最適化: f(x) + k · g(x) で c_0, c_1 を調整
// g(x) = c_d · x - m (有理側多項式)
// k を走査して polyScore を最小化
static void rotateOptimize(Poly& f, int d, const Int& m) {
    // g(x) = c_d · x - m → f(x) + k·g(x) で:
    //   c_0 → c_0 - k·m
    //   c_1 → c_1 + k·c_d
    //   (上位係数は不変)
    Int cd = f[d];

    // k の探索範囲: |c_0 - k·m| を最小化する k₀ の近傍
    // k₀ ≈ c_0 / m
    Int k0(0);
    if (!m.isZero() && m.bitLength() > 0) {
        k0 = f[0] / m;
    }

    double best_score = polyScore(f, d);
    Int best_k(0);

    // k₀ ± search_range で探索
    int search_range = 1000;
    for (int dk = -search_range; dk <= search_range; dk++) {
        Int k = k0 + Int(dk);

        // 仮の係数を計算
        Poly f2 = f;
        f2[0] = f[0] - k * m;
        f2[1] = f[1] + k * cd;

        double s = polyScore(f2, d);
        if (s < best_score) {
            best_score = s;
            best_k = k;
        }
    }

    // 最良の k を適用
    if (best_k != Int(0)) {
        f[0] = f[0] - best_k * m;
        f[1] = f[1] + best_k * cd;
    }
}

// 二次回転最適化: f(x) + (k₀ + k₁·x) · g(x)
// c_0 → c_0 - k₀·m, c_1 → c_1 + k₀·c_d - k₁·m, c_2 → c_2 + k₁·c_d
static void rotateOptimize2(Poly& f, int d, const Int& m) {
    if (d < 3) {
        rotateOptimize(f, d, m);
        return;
    }

    Int cd = f[d];
    // まず k₀ で c_0 を最小化
    Int k0_center(0);
    if (!m.isZero()) k0_center = f[0] / m;

    double best_score = polyScore(f, d);
    Int best_k0(0), best_k1(0);

    // 粗い探索: k₀ ± 500
    for (int dk0 = -500; dk0 <= 500; dk0 += 5) {
        Int k0 = k0_center + Int(dk0);
        Int new_c0 = f[0] - k0 * m;
        Int new_c1 = f[1] + k0 * cd;

        // k₁ で c_1 を最小化: k₁ ≈ new_c1 / m (ただし c_2 との
        // トレードオフ)
        Int k1_center(0);
        if (!m.isZero() && m.bitLength() > 1) k1_center = new_c1 / m;

        for (int dk1 = -20; dk1 <= 20; dk1++) {
            Int k1 = k1_center + Int(dk1);
            Poly f2 = f;
            f2[0] = new_c0;
            f2[1] = new_c1 - k1 * m;
            f2[2] = f[2] + k1 * cd;

            double s = polyScore(f2, d);
            if (s < best_score) {
                best_score = s;
                best_k0 = k0;
                best_k1 = k1;
            }
        }
    }

    // 精密探索: best_k0 ± 10
    Int k0_fine = best_k0;
    Int k1_fine = best_k1;
    for (int dk0 = -10; dk0 <= 10; dk0++) {
        Int k0 = k0_fine + Int(dk0);
        Int new_c0 = f[0] - k0 * m;
        Int new_c1 = f[1] + k0 * cd;

        Int k1c = k1_fine;
        for (int dk1 = -10; dk1 <= 10; dk1++) {
            Int k1 = k1c + Int(dk1);
            Poly f2 = f;
            f2[0] = new_c0;
            f2[1] = new_c1 - k1 * m;
            f2[2] = f[2] + k1 * cd;

            double s = polyScore(f2, d);
            if (s < best_score) {
                best_score = s;
                best_k0 = k0;
                best_k1 = k1;
            }
        }
    }

    // 最良を適用
    if (best_k0 != Int(0) || best_k1 != Int(0)) {
        f[0] = f[0] - best_k0 * m;
        f[1] = f[1] + best_k0 * cd - best_k1 * m;
        f[2] = f[2] + best_k1 * cd;
    }
}

// Kleinjung-style polynomial selection
// Smooth c_d 列挙 → 格子ベース m 探索 → 回転最適化
GNFSPoly kleinjungSelect(const Int& n, int d) {
    auto ipow = [](const Int& base, int exp) -> Int {
        Int r(1); for (int i = 0; i < exp; i++) r *= base; return r;
    };

    // c_d 上限: 次数に応じて調整
    uint64_t cd_max;
    switch (d) {
        case 3: cd_max = 10000; break;
        case 4: cd_max = 5000;  break;
        case 5: cd_max = 2000;  break;
        default: cd_max = 1000; break;
    }
    auto smooth_cds = enumerateSmoothCd(d, cd_max);

    // 射影根スコアでソート (高い順)
    std::vector<std::pair<double, uint64_t>> scored_cds;
    for (uint64_t cd : smooth_cds) {
        if (cd == 0) continue;
        scored_cds.push_back({projAlphaScore(cd), cd});
    }
    std::sort(scored_cds.begin(), scored_cds.end(),
              [](auto& a, auto& b) { return a.first > b.first; });

    // 上位候補に絞る (全探索は重いので)
    size_t max_cds = std::min(scored_cds.size(), (size_t)200);

    GNFSPoly best;
    best.d = d;
    best.score = 1e100;
    bool found = false;

    for (size_t ci = 0; ci < max_cds; ci++) {
        uint64_t cd_val = scored_cds[ci].second;
        Int cd(cd_val);

        // m₀ = floor((n / c_d)^{1/d})
        Int n_div = n / cd;
        if (n_div <= 1) continue;

        double nd = std::log2(std::max(1.0, n_div.toDouble())) / d;
        if (nd > 62.0) {
            // toDouble overflow 回避: Int の nthRoot を使う
            // 簡易版: 二分探索で m₀ を求める
            Int lo(1), hi = n_div;
            // hi を絞る: 2^(nd) の近似
            int bit_est = (int)(n_div.bitLength() / d);
            if (bit_est > 1) {
                hi = Int(1) << (bit_est + 2);
                lo = Int(1) << (bit_est - 1);
            }
            Int m0 = lo;
            for (int iter = 0; iter < 200; iter++) {
                Int mid = (lo + hi) / 2;
                if (mid == lo) break;
                Int md = ipow(mid, d);
                if (md <= n_div) { lo = mid; m0 = mid; }
                else hi = mid;
            }
            // m₀ の微調整
            while (ipow(m0 + 1, d) <= n_div) m0 += 1;

            if (m0 < 2) continue;

            // base-m 展開
            Poly f2(d + 1);
            Int rem = n;
            f2[d] = cd;
            rem -= cd * ipow(m0, d);
            if (rem.isNegative()) {
                m0 -= 1;
                rem = n - cd * ipow(m0, d);
            }
            if (rem.isNegative()) continue;

            for (int i = d - 1; i >= 1; i--) {
                Int mi = ipow(m0, i);
                if (mi.isZero()) break;
                f2[i] = rem / mi;
                rem -= f2[i] * mi;
            }
            f2[0] = rem;

            // 検証: f(m) = n
            {
                Int check(0), mp(1);
                for (int i = 0; i <= d; i++) { check += f2[i] * mp; mp *= m0; }
                if (check != n) continue;
            }

            // 回転最適化
            rotateOptimize2(f2, d, m0);

            // f(m) = n の再検証 (回転後)
            {
                Int check(0), mp(1);
                for (int i = 0; i <= d; i++) { check += f2[i] * mp; mp *= m0; }
                if (check != n) continue;
            }

            double s = polyScore(f2, d);
            if (s < best.score) {
                best.f = f2;
                best.m = m0;
                best.score = s;
                best.skewness = polySkewness(f2, d);
                found = true;
            }
        } else {
            Int m0 = Int(static_cast<uint64_t>(std::pow(2.0, nd)));
            if (m0 < 2) continue;

            // m^d · c_d ≤ n の保証
            for (int iter = 0; iter < 100; iter++) {
                Int md = ipow(m0, d) * cd;
                if (md > n) { m0 -= 1; continue; }
                Int md1 = ipow(m0 + 1, d) * cd;
                if (md1 <= n) { m0 += 1; continue; }
                break;
            }
            if (m0 < 2) continue;

            // ── 格子ベース探索: c_{d-1} を最小化 ──
            // base-m 展開で初期 c_{d-1} を取得
            Poly f_base(d + 1);
            {
                Int rem = n;
                f_base[d] = cd;
                rem -= cd * ipow(m0, d);
                if (rem.isNegative()) continue;
                for (int i = d - 1; i >= 1; i--) {
                    Int mi = ipow(m0, i);
                    if (mi.isZero()) break;
                    f_base[i] = rem / mi;
                    rem -= f_base[i] * mi;
                }
                f_base[0] = rem;

                // 検証
                Int check(0), mp(1);
                for (int i = 0; i <= d; i++) { check += f_base[i] * mp; mp *= m0; }
                if (check != n) continue;
            }

            // 格子 L: δm の変化で c_{d-1} がどう変わるか
            // f(m+δ) = n → c_{d-1} は c_{d-1} - cd·d·δ·m^{d-2} 近似
            // 格子基底: (1, -cd·d·m^{d-2}), (0, m^{d-1})
            // これをサイズ縮約して短ベクトル (δm, δc_{d-1}) を得る
            Int md2 = (d >= 2) ? ipow(m0, d - 2) : 1;
            Int md1 = ipow(m0, d - 1);

            Int e1x = 1;
            Int e1y = -(cd * d * md2);
            Int e2x = 0;
            Int e2y = md1;

            latticeReduce2D(e1x, e1y, e2x, e2y);

            // 短ベクトル e1 方向に ±K 個の候補を列挙
            int K = 200;
            for (int dk = -K; dk <= K; dk++) {
                Int delta_m = e1x * dk;
                Int m2 = m0 + delta_m;
                if (m2 < 2) continue;

                // base-m 展開で多項式を構成
                Poly f2(d + 1);
                f2[d] = cd;
                Int rem = n - cd * ipow(m2, d);
                if (rem.isNegative()) continue;

                bool valid = true;
                for (int i = d - 1; i >= 1; i--) {
                    Int mi = ipow(m2, i);
                    if (mi.isZero()) { valid = false; break; }
                    f2[i] = rem / mi;
                    rem -= f2[i] * mi;
                }
                if (!valid) continue;
                f2[0] = rem;

                // 検証: f(m2) = n
                {
                    Int check(0), mp(1);
                    for (int i = 0; i <= d; i++) { check += f2[i] * mp; mp *= m2; }
                    if (check != n) continue;
                }

                // 回転最適化
                rotateOptimize2(f2, d, m2);

                // 再検証
                {
                    Int check(0), mp(1);
                    for (int i = 0; i <= d; i++) { check += f2[i] * mp; mp *= m2; }
                    if (check != n) continue;
                }

                double s = polyScore(f2, d);
                if (s < best.score) {
                    best.f = f2;
                    best.m = m2;
                    best.score = s;
                    best.skewness = polySkewness(f2, d);
                    found = true;
                }
            }

            // e2 方向も少し探索
            for (int dk = -50; dk <= 50; dk++) {
                if (dk == 0) continue;
                Int delta_m = e2x * Int(dk);
                if (delta_m.isZero()) continue;  // e2x=0 なら m は変わらない
                Int m2 = m0 + delta_m;
                if (m2 < 2) continue;

                Poly f2(d + 1);
                f2[d] = cd;
                Int rem = n - cd * ipow(m2, d);
                if (rem.isNegative()) continue;

                bool valid = true;
                for (int i = d - 1; i >= 1; i--) {
                    Int mi = ipow(m2, i);
                    if (mi.isZero()) { valid = false; break; }
                    f2[i] = rem / mi;
                    rem -= f2[i] * mi;
                }
                if (!valid) continue;
                f2[0] = rem;

                Int check(0), mp(1);
                for (int i = 0; i <= d; i++) { check += f2[i] * mp; mp *= m2; }
                if (check != n) continue;

                rotateOptimize2(f2, d, m2);

                {
                    Int check2(0), mp2(1);
                    for (int i = 0; i <= d; i++) { check2 += f2[i] * mp2; mp2 *= m2; }
                    if (check2 != n) continue;
                }

                double s = polyScore(f2, d);
                if (s < best.score) {
                    best.f = f2;
                    best.m = m2;
                    best.score = s;
                    best.skewness = polySkewness(f2, d);
                    found = true;
                }
            }
        }
    }

    // フォールバック: 見つからなければ baseMSelect
    if (!found) {
        return baseMSelect(n, d);
    }

    return best;
}

// ============================================================
// Algebraic norm: N(a + bα) = (-b)^d · f(-a/b) (homogeneous)
// ============================================================

Int algNorm(int64_t a, int64_t b, const Poly& f, int d) {
    Int result(0);
    Int neg_a(-a), B(b);
    Int na_pow(1);
    std::vector<Int> bp(d + 1);
    bp[0] = 1;
    for (int i = 1; i <= d; i++) bp[i] = bp[i-1] * B;

    for (int i = 0; i <= d; i++) {
        result += f[i] * na_pow * bp[d - i];
        na_pow *= neg_a;
    }
    if (d % 2 == 1) result = -result;
    return result;
}

// f(r) mod p (uint64)
uint64_t evalFmodP(const Poly& f, int d, uint64_t r, uint64_t p) {
    uint64_t val = 0, rp = 1;
    for (int i = 0; i <= d; i++) {
        Int ci = f[i] % Int(p);
        if (ci.isNegative()) ci += Int(p);
        val = (val + ci.toUInt64() * rp) % p;
        rp = rp * r % p;
    }
    return val;
}

// ============================================================
// Factor bases
// ============================================================

struct AlgFBEntry { uint64_t p, r; };

std::vector<uint64_t> buildRatFB(uint64_t bound) {
    std::vector<uint64_t> fb = {2};
    for (uint64_t p = 3; p <= bound; p += 2) {
        bool ok = true;
        for (uint64_t d = 3; d*d <= p; d += 2) if (p % d == 0) { ok = false; break; }
        if (ok) fb.push_back(p);
    }
    return fb;
}

// ============================================================
// GF(p) 上の多項式演算 (Cantor-Zassenhaus 用)
// 係数は uint64_t, p < 2^32
// ============================================================

using PolyGFp = std::vector<uint64_t>;

// 多項式の次数 (零多項式は -1)
static int polyDeg(const PolyGFp& a) {
    int d = (int)a.size() - 1;
    while (d >= 0 && a[d] == 0) d--;
    return d;
}

// 正規化 (上位の 0 係数を除去)
static PolyGFp polyTrim(PolyGFp a) {
    while (a.size() > 1 && a.back() == 0) a.pop_back();
    return a;
}

// 多項式乗算 mod p
static PolyGFp polyMulGFp(const PolyGFp& a, const PolyGFp& b, uint64_t p) {
    int da = polyDeg(a), db = polyDeg(b);
    if (da < 0 || db < 0) return {0};
    PolyGFp r(da + db + 1, 0);
    for (int i = 0; i <= da; i++)
        for (int j = 0; j <= db; j++)
            r[i + j] = (r[i + j] + a[i] * b[j]) % p;
    return polyTrim(r);
}

// 多項式剰余 a mod b mod p
static PolyGFp polyModGFp(PolyGFp a, const PolyGFp& b, uint64_t p) {
    int db = polyDeg(b);
    if (db < 0) return {0};
    uint64_t lc_inv = gnfs_modinv(b[db], p);
    while (true) {
        int da = polyDeg(a);
        if (da < db) break;
        uint64_t scale = a[da] * lc_inv % p;
        for (int i = 0; i <= db; i++)
            a[da - db + i] = (a[da - db + i] + p - scale * b[i] % p) % p;
    }
    return polyTrim(a);
}

// 多項式 GCD mod p
static PolyGFp polyGcdGFp(PolyGFp a, PolyGFp b, uint64_t p) {
    while (polyDeg(b) >= 0) {
        a = polyModGFp(a, b, p);
        std::swap(a, b);
    }
    if (polyDeg(a) >= 0) {
        uint64_t lc_inv = gnfs_modinv(a[polyDeg(a)], p);
        for (auto& c : a) c = c * lc_inv % p;
    }
    return polyTrim(a);
}

// x^e mod g mod p (repeated squaring)
static PolyGFp polyPowModGFp(PolyGFp base, uint64_t e, const PolyGFp& g, uint64_t p) {
    PolyGFp result = {1};
    base = polyModGFp(base, g, p);
    while (e > 0) {
        if (e & 1) result = polyModGFp(polyMulGFp(result, base, p), g, p);
        base = polyModGFp(polyMulGFp(base, base, p), g, p);
        e >>= 1;
    }
    return result;
}

// Cantor-Zassenhaus: f(x) mod p の全根を返す
// f は monic (正規化済み) を想定
static std::vector<uint64_t> findRootsModP(const Poly& f_int, int d, uint64_t p) {
    // f を GF(p) 多項式に変換
    PolyGFp fp(d + 1);
    for (int i = 0; i <= d; i++) {
        Int ci = f_int[i] % Int(p);
        if (ci.isNegative()) ci += Int(p);
        fp[i] = ci.toUInt64();
    }
    fp = polyTrim(fp);

    // p == 2 は直接評価
    if (p == 2) {
        std::vector<uint64_t> roots;
        for (uint64_t r = 0; r < 2; r++)
            if (evalFmodP(f_int, d, r, p) == 0) roots.push_back(r);
        return roots;
    }

    // Step 1: gcd(x^p - x, f) で根を持つ因子を抽出
    // x^p mod f
    PolyGFp xp = polyPowModGFp({0, 1}, p, fp, p);
    // x^p - x
    if (xp.size() < 2) xp.resize(2, 0);
    xp[1] = (xp[1] + p - 1) % p;
    PolyGFp g0 = polyGcdGFp(fp, xp, p);

    int nroots = polyDeg(g0);
    if (nroots <= 0) return {};

    // Step 2: 再帰的に 1 次因子に分解
    std::vector<PolyGFp> factors;
    std::vector<PolyGFp> work = {g0};
    std::mt19937_64 rng(p * 1000003ULL + 42);

    while (!work.empty()) {
        PolyGFp h = work.back(); work.pop_back();
        int dh = polyDeg(h);
        if (dh == 1) {
            factors.push_back(h);
            continue;
        }
        if (dh == 0) continue;

        // ランダム分割: gcd((a + x)^{(p-1)/2} - 1, h)
        for (int attempt = 0; attempt < 100; attempt++) {
            uint64_t a = rng() % p;
            PolyGFp ax = {a, 1}; // a + x
            PolyGFp t = polyPowModGFp(ax, (p - 1) / 2, h, p);
            // t - 1
            t[0] = (t[0] + p - 1) % p;
            PolyGFp g = polyGcdGFp(t, h, p);
            int dg = polyDeg(g);
            if (dg > 0 && dg < dh) {
                work.push_back(g);
                work.push_back(polyModGFp(h, g, p) /* h/g */);
                // h / g
                PolyGFp& last = work.back();
                last = {0}; // 再計算
                // h を g で割る
                PolyGFp quot(dh - dg + 1, 0);
                PolyGFp rem = h;
                uint64_t glc_inv = gnfs_modinv(g[dg], p);
                for (int i = dh; i >= dg; i--) {
                    int dr = polyDeg(rem);
                    if (dr < dg) break;
                    uint64_t q = rem[dr] * glc_inv % p;
                    quot[dr - dg] = q;
                    for (int j = 0; j <= dg; j++)
                        rem[dr - dg + j] = (rem[dr - dg + j] + p - q * g[j] % p) % p;
                }
                last = polyTrim(quot);
                break;
            }
        }
    }

    // 1 次因子 (x - r) から根を抽出
    std::vector<uint64_t> roots;
    for (auto& fac : factors) {
        if (polyDeg(fac) != 1) continue;
        // x + c → 根は -c mod p
        uint64_t r = (p - fac[0]) % p;
        roots.push_back(r);
    }
    std::sort(roots.begin(), roots.end());
    return roots;
}

std::vector<AlgFBEntry> buildAlgFB(const Poly& f, int d, const std::vector<uint64_t>& rfb) {
    std::vector<AlgFBEntry> afb;
    for (uint64_t p : rfb) {
        auto roots = findRootsModP(f, d, p);
        for (uint64_t r : roots)
            afb.push_back({p, r});
    }
    return afb;
}

// ============================================================
// Sieve + relation collection
// ============================================================

struct GNFSRel {
    int64_t a, b;
    std::vector<int> rexp, aexp;
    bool rat_neg;
};

// Large prime partial relation
// カテゴリ:
//   RF: rat_lp_count=1, alg_lp=0 (rational 1-LP, algebraic smooth)
//   AF: rat_lp_count=0, alg_lp>0 (rational smooth, algebraic 1-LP)
//   2RF: rat_lp_count=2, alg_lp=0 (rational 2-LP, algebraic smooth)
struct GNFSPartial {
    int64_t a, b;
    std::vector<int> rexp, aexp;
    bool rat_neg;
    uint64_t rat_lp[2];     // rational large primes (0 = none)
    int rat_lp_count;        // 0, 1, or 2
    uint64_t alg_lp;         // algebraic large prime (0 = none)
};

std::vector<GNFSRel> gnfsCollect(
    const GNFSPoly& poly, const std::vector<uint64_t>& rfb,
    const std::vector<AlgFBEntry>& afb, const GNFSParams& par)
{
    int rn = (int)rfb.size(), an = (int)afb.size();
    int target = rn + an + 10;
    std::vector<GNFSRel> rels;
    std::vector<GNFSPartial> partials;
    uint64_t lp_bound = rfb.back() * rfb.back();       // 1-LP 上限
    uint64_t lp2_bound = lp_bound;                      // 2-LP: 各素因数の上限
    // 2-LP 余因子上限: lp2_bound² (両方 ≤ lp2_bound)
    // overflow 回避: lp2_bound ≤ 2^32 なら lp2_bound² ≤ 2^64
    uint64_t lp2_cofactor_bound = (lp2_bound <= 0xFFFFFFFFULL)
        ? lp2_bound * lp2_bound : 0xFFFFFFFFFFFFFFFFULL;
    int A = par.sieve_a, slen = 2*A + 1;

    // m mod p for rational sieve
    std::vector<uint64_t> mmodp(rn);
    for (int i = 0; i < rn; i++) {
        Int mm = poly.m % Int(rfb[i]);
        if (mm.isNegative()) mm += Int(rfb[i]);
        mmodp[i] = mm.toUInt64();
    }

    std::vector<float> rlog(rn), alog(an);
    for (int i = 0; i < rn; i++) rlog[i] = (float)std::log2((double)rfb[i]);
    for (int i = 0; i < an; i++) alog[i] = (float)std::log2((double)afb[i].p);

    constexpr int GNFS_BLOCK = 8192;
    std::vector<float> sr(GNFS_BLOCK), sa(GNFS_BLOCK);

    for (int64_t b = 1; b <= par.sieve_b && (int)rels.size() < target; b++) {
        // 篩開始位置を事前計算
        std::vector<int64_t> rstart(rn), astart(an);
        for (int i = 0; i < rn; i++) {
            uint64_t p = rfb[i];
            uint64_t bmodp = (uint64_t)(b % (int64_t)p);
            uint64_t start = (p - bmodp * mmodp[i] % p) % p;
            rstart[i] = (int64_t)((start + (uint64_t)(A % (int64_t)p)) % p);
        }
        for (int i = 0; i < an; i++) {
            uint64_t p = afb[i].p, r = afb[i].r;
            uint64_t bmodp = (uint64_t)(b % (int64_t)p);
            uint64_t start = (p - bmodp * r % p) % p;
            astart[i] = (int64_t)((start + (uint64_t)(A % (int64_t)p)) % p);
        }

        // Threshold
        double bm = std::abs((double)b * poly.m.toDouble());
        float rthr = (float)(std::log2(std::max(bm, 1.0) + A) - par.thresh_adj);
        float athr = (float)((double)poly.d * std::log2(std::max((double)A, 1.0)) - par.thresh_adj);
        if (rthr < 3) rthr = 3;
        if (athr < 3) athr = 3;

        // ブロック篩
        for (int blk_start = 0; blk_start < slen && (int)rels.size() < target;
             blk_start += GNFS_BLOCK) {
            int blk_end = std::min(blk_start + GNFS_BLOCK, slen);
            int blk_len = blk_end - blk_start;
            std::fill(sr.begin(), sr.begin() + blk_len, 0.0f);
            std::fill(sa.begin(), sa.begin() + blk_len, 0.0f);

            // Rational sieve: a ≡ -b·m (mod p)
            for (int i = 0; i < rn; i++) {
                uint64_t p = rfb[i];
                int64_t j = rstart[i];
                if (j < blk_start) {
                    int64_t skip = (blk_start - j + (int64_t)p - 1) / (int64_t)p;
                    j += skip * (int64_t)p;
                }
                for (; j < blk_end; j += (int64_t)p)
                    sr[j - blk_start] += rlog[i];
            }

            // Algebraic sieve: a ≡ -b·r (mod p)
            for (int i = 0; i < an; i++) {
                uint64_t p = afb[i].p;
                int64_t j = astart[i];
                if (j < blk_start) {
                    int64_t skip = (blk_start - j + (int64_t)p - 1) / (int64_t)p;
                    j += skip * (int64_t)p;
                }
                for (; j < blk_end; j += (int64_t)p)
                    sa[j - blk_start] += alog[i];
            }

            for (int idx_local = 0; idx_local < blk_len && (int)rels.size() < target; idx_local++) {
                if (sr[idx_local] < rthr || sa[idx_local] < athr) continue;
                int idx = blk_start + idx_local;
            int64_t a = idx - A;
            if (a == 0) continue;
            if (std::gcd(std::abs(a), b) != 1) continue;

            // Trial divide rational norm |a + b·m|
            Int rv = Int(a) + Int(b) * poly.m;
            bool rneg = rv.isNegative();
            Int ra = rneg ? -rv : rv;
            std::vector<int> re(rn, 0);
            bool r_smooth;
            uint64_t r_rem64 = 0;
            if (ra.bitLength() <= 64) {
                r_rem64 = ra.toUInt64();
                for (int i = 0; i < rn && r_rem64 > 1; i++) {
                    uint64_t p = rfb[i];
                    while (r_rem64 % p == 0) { r_rem64 /= p; re[i]++; }
                }
                r_smooth = (r_rem64 == 1);
            } else {
                Int rrem = ra;
                for (int i = 0; i < rn; i++) {
                    Int p(rfb[i]);
                    while ((rrem % p).isZero()) { rrem /= p; re[i]++; }
                }
                r_smooth = (rrem == 1);
                if (!r_smooth && rrem.bitLength() <= 64) r_rem64 = rrem.toUInt64();
            }
            // Rational cofactor 分類
            int r_lp_count = 0;
            uint64_t r_lps[2] = {0, 0};
            if (!r_smooth) {
                if (r_rem64 <= 1) continue;
                if (r_rem64 <= lp_bound) {
                    // 1-LP 候補 (素数判定は省略: FB 以下の因子なし → 高確率で素数)
                    r_lp_count = 1;
                    r_lps[0] = r_rem64;
                } else if (r_rem64 <= lp2_cofactor_bound) {
                    // 2-LP 候補: 余因子が 2 つの素数の積か判定
                    if (isPrime64(r_rem64)) continue;  // 単一の大素数 → 棄却
                    uint64_t f1 = pollardRho64(r_rem64);
                    if (f1 == 0) continue;
                    uint64_t f2 = r_rem64 / f1;
                    if (f1 > f2) std::swap(f1, f2);
                    if (f2 > lp2_bound) continue;  // 片方が大きすぎる
                    r_lp_count = 2;
                    r_lps[0] = f1; r_lps[1] = f2;
                } else {
                    continue;  // 余因子が大きすぎる
                }
            }

            // Trial divide algebraic norm
            Int norm = algNorm(a, b, poly.f, poly.d);
            if (norm.isZero()) continue;
            Int na = norm.isNegative() ? -norm : norm;
            std::vector<int> ae(an, 0);
            bool a_smooth;
            uint64_t a_rem64 = 0;
            if (na.bitLength() <= 64) {
                a_rem64 = na.toUInt64();
                for (int i = 0; i < an && a_rem64 > 1; i++) {
                    uint64_t p = afb[i].p;
                    while (a_rem64 % p == 0) { a_rem64 /= p; ae[i]++; }
                }
                a_smooth = (a_rem64 == 1);
            } else {
                Int arem = na;
                for (int i = 0; i < an; i++) {
                    Int p(afb[i].p);
                    while ((arem % p).isZero()) { arem /= p; ae[i]++; }
                }
                a_smooth = (arem == 1);
                if (!a_smooth && arem.bitLength() <= 64) a_rem64 = arem.toUInt64();
            }

            // Algebraic cofactor 分類
            uint64_t a_lp = 0;
            if (!a_smooth) {
                // Algebraic 1-LP: 余因子が LP 範囲内なら受理
                // rat 側が既に LP を持つ場合 (r_lp_count > 0) は棄却
                // (rat 1-LP + alg 1-LP の組は実用的にペアリング困難)
                if (r_lp_count > 0) continue;
                if (a_rem64 <= 1 || a_rem64 > lp_bound) continue;
                a_lp = a_rem64;
            }

            if (r_lp_count == 0 && a_lp == 0) {
                // FF: 両側完全 smooth
                rels.push_back({a, b, re, ae, rneg});
            } else {
                // Partial relation
                GNFSPartial part;
                part.a = a; part.b = b;
                part.rexp = re; part.aexp = ae;
                part.rat_neg = rneg;
                part.rat_lp[0] = r_lps[0]; part.rat_lp[1] = r_lps[1];
                part.rat_lp_count = r_lp_count;
                part.alg_lp = a_lp;
                partials.push_back(std::move(part));
            }
            }  // smooth 候補ループ end
        }  // ブロックループ end
    }

    // ── Phase 1: Rational 1-LP ペアリング (RF 型) ──
    {
        std::unordered_map<uint64_t, size_t> lp_map;
        for (size_t i = 0; i < partials.size(); i++) {
            if (partials[i].rat_lp_count != 1 || partials[i].alg_lp != 0) continue;
            auto it = lp_map.find(partials[i].rat_lp[0]);
            if (it != lp_map.end()) {
                auto& p1 = partials[it->second];
                auto& p2 = partials[i];
                GNFSRel combined;
                combined.a = p1.a; combined.b = p1.b;
                combined.rexp.resize(rn, 0);
                combined.aexp.resize(an, 0);
                for (int j = 0; j < rn; j++)
                    combined.rexp[j] = p1.rexp[j] + p2.rexp[j];
                for (int j = 0; j < an; j++)
                    combined.aexp[j] = p1.aexp[j] + p2.aexp[j];
                combined.rat_neg = (p1.rat_neg != p2.rat_neg);
                rels.push_back(combined);
            } else {
                lp_map[partials[i].rat_lp[0]] = i;
            }
        }
    }

    // ── Phase 2: Algebraic 1-LP ペアリング (AF 型) ──
    {
        std::unordered_map<uint64_t, size_t> alp_map;
        for (size_t i = 0; i < partials.size(); i++) {
            if (partials[i].rat_lp_count != 0 || partials[i].alg_lp == 0) continue;
            auto it = alp_map.find(partials[i].alg_lp);
            if (it != alp_map.end()) {
                auto& p1 = partials[it->second];
                auto& p2 = partials[i];
                GNFSRel combined;
                combined.a = p1.a; combined.b = p1.b;
                combined.rexp.resize(rn, 0);
                combined.aexp.resize(an, 0);
                for (int j = 0; j < rn; j++)
                    combined.rexp[j] = p1.rexp[j] + p2.rexp[j];
                for (int j = 0; j < an; j++)
                    combined.aexp[j] = p1.aexp[j] + p2.aexp[j];
                combined.rat_neg = (p1.rat_neg != p2.rat_neg);
                rels.push_back(combined);
            } else {
                alp_map[partials[i].alg_lp] = i;
            }
        }
    }

    // ── Phase 3: Rational 2-LP グラフサイクル合成 (2RF 型) ──
    // 頂点 = large prime, 辺 = 2-LP partial relation
    // サイクル内の全 partial を XOR すると LP が偶数回現れ GF(2) で消える
    {
        // 2RF partial のインデックスを収集
        std::vector<size_t> idx_2lp;
        for (size_t i = 0; i < partials.size(); i++) {
            if (partials[i].rat_lp_count == 2 && partials[i].alg_lp == 0)
                idx_2lp.push_back(i);
        }

        if (!idx_2lp.empty()) {
            // 隣接リスト構築: node → [(neighbor, edge_index_in_idx_2lp)]
            std::unordered_map<uint64_t, std::vector<std::pair<uint64_t, size_t>>> adj;
            for (size_t ei = 0; ei < idx_2lp.size(); ei++) {
                auto& p = partials[idx_2lp[ei]];
                uint64_t u = p.rat_lp[0], v = p.rat_lp[1];
                if (u == v) {
                    // 完全平方余因子 → LP は偶数指数 → そのまま full relation
                    GNFSRel combined;
                    combined.a = p.a; combined.b = p.b;
                    combined.rexp = p.rexp; combined.aexp = p.aexp;
                    combined.rat_neg = p.rat_neg;
                    rels.push_back(combined);
                    continue;
                }
                adj[u].push_back({v, ei});
                adj[v].push_back({u, ei});
            }

            // BFS で全域森を構築
            std::unordered_map<uint64_t, uint64_t> parent;
            std::unordered_map<uint64_t, size_t> parent_edge;
            std::unordered_map<uint64_t, int> depth;
            std::unordered_set<size_t> tree_edges;

            for (auto& [start, _] : adj) {
                if (parent.count(start)) continue;
                parent[start] = start;
                depth[start] = 0;
                std::queue<uint64_t> q;
                q.push(start);
                while (!q.empty()) {
                    uint64_t u = q.front(); q.pop();
                    for (auto& [v, ei] : adj[u]) {
                        if (!parent.count(v)) {
                            parent[v] = u;
                            parent_edge[v] = ei;
                            depth[v] = depth[u] + 1;
                            tree_edges.insert(ei);
                            q.push(v);
                        }
                    }
                }
            }

            // 非木辺 → サイクル検出 & 合成
            for (size_t ei = 0; ei < idx_2lp.size(); ei++) {
                if (tree_edges.count(ei)) continue;
                auto& pe = partials[idx_2lp[ei]];
                uint64_t u = pe.rat_lp[0], v = pe.rat_lp[1];
                if (u == v) continue;  // 既処理

                // u, v から LCA まで木をたどり、パス上の辺を収集
                std::vector<size_t> cycle_edges;
                cycle_edges.push_back(ei);  // 非木辺自身

                uint64_t a = u, b = v;
                // 深さを揃える
                while (depth[a] > depth[b]) {
                    cycle_edges.push_back(parent_edge[a]);
                    a = parent[a];
                }
                while (depth[b] > depth[a]) {
                    cycle_edges.push_back(parent_edge[b]);
                    b = parent[b];
                }
                // LCA まで同時にたどる
                while (a != b) {
                    cycle_edges.push_back(parent_edge[a]);
                    cycle_edges.push_back(parent_edge[b]);
                    a = parent[a];
                    b = parent[b];
                }

                // サイクル内の全 partial を XOR 合成
                GNFSRel combined;
                combined.a = partials[idx_2lp[cycle_edges[0]]].a;
                combined.b = partials[idx_2lp[cycle_edges[0]]].b;
                combined.rexp.assign(rn, 0);
                combined.aexp.assign(an, 0);
                combined.rat_neg = false;
                for (size_t ce : cycle_edges) {
                    auto& cp = partials[idx_2lp[ce]];
                    for (int j = 0; j < rn; j++) combined.rexp[j] += cp.rexp[j];
                    for (int j = 0; j < an; j++) combined.aexp[j] += cp.aexp[j];
                    combined.rat_neg = (combined.rat_neg != cp.rat_neg);
                }
                rels.push_back(combined);
            }
        }
    }

    return rels;
}

// ============================================================
// Lattice sieve (special-q)
// ============================================================
// 大きな素数イデアル (q, s) を "special-q" として選択し、
// 格子 L_q = {(a,b) : a + b·s ≡ 0 mod q} に制限して篩う。
// q は代数側で自動的に smooth (1 個分の因子を「無料」で獲得)。
// 2D Lagrange 格子基底縮約で短い基底ベクトルを計算。

void gnfsLatticeSieve(
    const GNFSPoly& poly, const std::vector<uint64_t>& rfb,
    const std::vector<AlgFBEntry>& afb, const GNFSParams& par,
    std::vector<GNFSRel>& rels, std::vector<GNFSPartial>& partials,
    uint64_t lp_bound, int target)
{
    int rn = (int)rfb.size(), an = (int)afb.size();

    std::vector<float> rlog(rn), alog(an);
    for (int i = 0; i < rn; i++) rlog[i] = (float)std::log2((double)rfb[i]);
    for (int i = 0; i < an; i++) alog[i] = (float)std::log2((double)afb[i].p);

    std::vector<uint64_t> mmodp(rn);
    for (int i = 0; i < rn; i++) {
        Int mm = poly.m % Int(rfb[i]);
        if (mm.isNegative()) mm += Int(rfb[i]);
        mmodp[i] = mm.toUInt64();
    }

    // special-q 範囲: FB 上限の 1-4 倍
    uint64_t q_min = rfb.back() + 2;
    uint64_t q_max = std::min(q_min * 4, q_min + 50000ULL);

    for (uint64_t q = q_min | 1; q <= q_max && (int)rels.size() < target; q += 2) {
        if (!isPrime64(q)) continue;
        auto qroots = findRootsModP(poly.f, poly.d, q);
        if (qroots.empty()) continue;

        for (uint64_t s : qroots) {
            if ((int)rels.size() >= target) break;

            // 格子基底: {(q, 0), (-(int64_t)s, 1)}
            int64_t e1x = (int64_t)q, e1y = 0;
            int64_t e2x = -(int64_t)s, e2y = 1;

            // 2D Lagrange 基底縮約
            auto n2 = [](int64_t x, int64_t y) -> double {
                return (double)x * x + (double)y * y;
            };
            for (int it = 0; it < 100; it++) {
                if (n2(e1x, e1y) > n2(e2x, e2y)) {
                    std::swap(e1x, e2x); std::swap(e1y, e2y);
                }
                double dot = (double)e1x * e2x + (double)e1y * e2y;
                double nn = n2(e1x, e1y);
                if (nn < 1.0) break;
                int64_t mu = (int64_t)std::round(dot / nn);
                if (mu == 0) break;
                e2x -= mu * e1x; e2y -= mu * e1y;
            }
            if (n2(e1x, e1y) > n2(e2x, e2y)) {
                std::swap(e1x, e2x); std::swap(e1y, e2y);
            }

            // 篩範囲: 格子座標 (u, v) ∈ [-U, U] × [1, V]
            double e1len = std::sqrt(n2(e1x, e1y));
            double e2len = std::sqrt(n2(e2x, e2y));
            int U = std::min((int)(par.sieve_a / std::max(e1len, 1.0)), 50000);
            int V = std::min((int)(par.sieve_b / std::max(e2len, 1.0)), 500);
            if (U < 10 || V < 1) continue;
            int slen = 2 * U + 1;

            // 各素数の格子内篩位置を事前計算
            // a + b·m ≡ 0 mod p → u·(e1x+m·e1y) + v·(e2x+m·e2y) ≡ 0 mod p
            // → u ≡ -v·C2·C1^{-1} mod p  (C1 = e1x+m·e1y mod p, C2 = e2x+m·e2y mod p)
            std::vector<uint64_t> rC1(rn), rC2(rn);
            std::vector<bool> rC1zero(rn, false);
            for (int i = 0; i < rn; i++) {
                uint64_t p = rfb[i];
                uint64_t ex = (uint64_t)((e1x % (int64_t)p + (int64_t)p) % (int64_t)p);
                uint64_t ey = (uint64_t)((e1y % (int64_t)p + (int64_t)p) % (int64_t)p);
                rC1[i] = (ex + mmodp[i] * ey) % p;
                ex = (uint64_t)((e2x % (int64_t)p + (int64_t)p) % (int64_t)p);
                ey = (uint64_t)((e2y % (int64_t)p + (int64_t)p) % (int64_t)p);
                rC2[i] = (ex + mmodp[i] * ey) % p;
                rC1zero[i] = (rC1[i] == 0);
            }

            // 代数側: a + b·r ≡ 0 mod p
            std::vector<uint64_t> aC1(an), aC2(an);
            std::vector<bool> aC1zero(an, false);
            for (int i = 0; i < an; i++) {
                uint64_t p = afb[i].p, r = afb[i].r;
                uint64_t ex = (uint64_t)((e1x % (int64_t)p + (int64_t)p) % (int64_t)p);
                uint64_t ey = (uint64_t)((e1y % (int64_t)p + (int64_t)p) % (int64_t)p);
                aC1[i] = (ex + r * ey % p) % p;
                ex = (uint64_t)((e2x % (int64_t)p + (int64_t)p) % (int64_t)p);
                ey = (uint64_t)((e2y % (int64_t)p + (int64_t)p) % (int64_t)p);
                aC2[i] = (ex + r * ey % p) % p;
                aC1zero[i] = (aC1[i] == 0);
            }

            constexpr int LSBLOCK = 8192;
            std::vector<float> sr(LSBLOCK), sa(LSBLOCK);

            for (int64_t v = 1; v <= V && (int)rels.size() < target; v++) {
                // ブロック篩 in u-direction
                for (int blk = -U; blk < U + 1 && (int)rels.size() < target; blk += LSBLOCK) {
                    int blk_end = std::min(blk + LSBLOCK, U + 1);
                    int blk_len = blk_end - blk;
                    std::fill(sr.begin(), sr.begin() + blk_len, 0.0f);
                    std::fill(sa.begin(), sa.begin() + blk_len, 0.0f);

                    // Rational sieve
                    for (int i = 0; i < rn; i++) {
                        uint64_t p = rfb[i];
                        if (rC1zero[i]) {
                            uint64_t vp = (uint64_t)((v % (int64_t)p + (int64_t)p) % (int64_t)p);
                            if ((rC2[i] * vp) % p == 0)
                                for (int j = 0; j < blk_len; j++) sr[j] += rlog[i];
                            continue;
                        }
                        uint64_t vp = (uint64_t)((v % (int64_t)p + (int64_t)p) % (int64_t)p);
                        uint64_t u0p = (p - rC2[i] * vp % p) % p * gnfs_modinv(rC1[i], p) % p;
                        // u0p は u mod p の篩位置 → blk 内で対応する位置へ
                        int64_t u_off = (int64_t)u0p - ((int64_t)blk % (int64_t)p + (int64_t)p) % (int64_t)p;
                        if (u_off < 0) u_off += (((-u_off) + (int64_t)p - 1) / (int64_t)p) * (int64_t)p;
                        for (int64_t j = u_off; j < blk_len; j += (int64_t)p)
                            sr[j] += rlog[i];
                    }

                    // Algebraic sieve
                    for (int i = 0; i < an; i++) {
                        uint64_t p = afb[i].p;
                        if (aC1zero[i]) {
                            uint64_t vp = (uint64_t)((v % (int64_t)p + (int64_t)p) % (int64_t)p);
                            if ((aC2[i] * vp) % p == 0)
                                for (int j = 0; j < blk_len; j++) sa[j] += alog[i];
                            continue;
                        }
                        uint64_t vp = (uint64_t)((v % (int64_t)p + (int64_t)p) % (int64_t)p);
                        uint64_t u0p = (p - aC2[i] * vp % p) % p * gnfs_modinv(aC1[i], p) % p;
                        int64_t u_off = (int64_t)u0p - ((int64_t)blk % (int64_t)p + (int64_t)p) % (int64_t)p;
                        if (u_off < 0) u_off += (((-u_off) + (int64_t)p - 1) / (int64_t)p) * (int64_t)p;
                        for (int64_t j = u_off; j < blk_len; j += (int64_t)p)
                            sa[j] += alog[i];
                    }

                    // special-q の寄与を加算 (代数側に log2(q) を追加)
                    float qlog = (float)std::log2((double)q);
                    for (int j = 0; j < blk_len; j++) sa[j] += qlog;

                    // Smooth 候補の検証
                    float rthr = (float)(std::log2(std::max((double)par.sieve_a, 1.0)) - par.thresh_adj + 5);
                    float athr = (float)((double)poly.d * std::log2(std::max((double)par.sieve_a, 1.0)) - par.thresh_adj + 5);
                    if (rthr < 3) rthr = 3;
                    if (athr < 3) athr = 3;

                    for (int j = 0; j < blk_len && (int)rels.size() < target; j++) {
                        if (sr[j] < rthr || sa[j] < athr) continue;
                        int64_t u = blk + j;
                        int64_t a = u * e1x + v * e2x;
                        int64_t b = u * e1y + v * e2y;
                        if (b <= 0 || a == 0) continue;
                        if (std::gcd(std::abs(a), b) != 1) continue;

                        // Rational trial division
                        Int rv = Int(a) + Int(b) * poly.m;
                        bool rneg = rv.isNegative();
                        Int ra = rneg ? -rv : rv;
                        std::vector<int> re(rn, 0);
                        bool r_smooth;
                        uint64_t r_rem64 = 0;
                        if (ra.bitLength() <= 64) {
                            r_rem64 = ra.toUInt64();
                            for (int i2 = 0; i2 < rn && r_rem64 > 1; i2++) {
                                uint64_t p = rfb[i2];
                                while (r_rem64 % p == 0) { r_rem64 /= p; re[i2]++; }
                            }
                            r_smooth = (r_rem64 == 1);
                        } else {
                            Int rrem = ra;
                            for (int i2 = 0; i2 < rn; i2++) {
                                Int p(rfb[i2]);
                                while ((rrem % p).isZero()) { rrem /= p; re[i2]++; }
                            }
                            r_smooth = (rrem == 1);
                            if (!r_smooth && rrem.bitLength() <= 64) r_rem64 = rrem.toUInt64();
                        }
                        if (!r_smooth && (r_rem64 <= 1 || r_rem64 > lp_bound)) continue;

                        // Algebraic trial division
                        Int norm = algNorm(a, b, poly.f, poly.d);
                        if (norm.isZero()) continue;
                        Int na = norm.isNegative() ? -norm : norm;
                        std::vector<int> ae(an, 0);
                        bool a_smooth;
                        uint64_t a_rem64 = 0;
                        if (na.bitLength() <= 64) {
                            a_rem64 = na.toUInt64();
                            // special-q の因子を除去
                            while (a_rem64 % q == 0) a_rem64 /= q;
                            for (int i2 = 0; i2 < an && a_rem64 > 1; i2++) {
                                uint64_t p = afb[i2].p;
                                while (a_rem64 % p == 0) { a_rem64 /= p; ae[i2]++; }
                            }
                            a_smooth = (a_rem64 == 1);
                        } else {
                            Int arem = na;
                            Int Q((uint64_t)q);
                            while ((arem % Q).isZero()) arem /= Q;
                            for (int i2 = 0; i2 < an; i2++) {
                                Int p(afb[i2].p);
                                while ((arem % p).isZero()) { arem /= p; ae[i2]++; }
                            }
                            a_smooth = (arem == 1);
                            if (!a_smooth && arem.bitLength() <= 64) a_rem64 = arem.toUInt64();
                        }
                        if (!a_smooth) continue;

                        if (r_smooth) {
                            rels.push_back({a, b, re, ae, rneg});
                        } else {
                            GNFSPartial part;
                            part.a = a; part.b = b;
                            part.rexp = re; part.aexp = ae;
                            part.rat_neg = rneg;
                            part.rat_lp[0] = r_rem64; part.rat_lp[1] = 0;
                            part.rat_lp_count = 1;
                            part.alg_lp = 0;
                            partials.push_back(std::move(part));
                        }
                    }
                }
            }
        }
    }
}

// ============================================================
// 関係のファイル入出力
// ============================================================

bool gnfsSaveRelations(const std::string& path,
    const Int& n, const GNFSPoly& poly,
    int rn, int an,
    const std::vector<GNFSRel>& rels)
{
    std::ofstream ofs(path);
    if (!ofs) return false;
    ofs << "GNFS_RELATIONS_V1\n";
    ofs << "n " << n.toString() << "\n";
    ofs << "poly " << poly.d;
    for (int i = 0; i <= poly.d; i++) ofs << " " << poly.f[i].toString();
    ofs << "\n";
    ofs << "m " << poly.m.toString() << "\n";
    ofs << "rfb_size " << rn << "\n";
    ofs << "afb_size " << an << "\n";
    ofs << "nrels " << rels.size() << "\n";
    for (auto& r : rels) {
        ofs << r.a << " " << r.b;
        for (int j = 0; j < rn; j++) ofs << " " << r.rexp[j];
        for (int j = 0; j < an; j++) ofs << " " << r.aexp[j];
        ofs << " " << (r.rat_neg ? 1 : 0) << "\n";
    }
    return ofs.good();
}

bool gnfsLoadRelations(const std::string& path,
    int rn, int an,
    std::vector<GNFSRel>& rels)
{
    std::ifstream ifs(path);
    if (!ifs) return false;
    std::string line;
    int nrels = 0;
    while (std::getline(ifs, line)) {
        if (line.substr(0, 5) == "nrels") {
            nrels = std::stoi(line.substr(6));
            break;
        }
    }
    rels.reserve(rels.size() + nrels);
    for (int i = 0; i < nrels && std::getline(ifs, line); i++) {
        std::istringstream iss(line);
        GNFSRel r;
        iss >> r.a >> r.b;
        r.rexp.resize(rn); r.aexp.resize(an);
        for (int j = 0; j < rn; j++) iss >> r.rexp[j];
        for (int j = 0; j < an; j++) iss >> r.aexp[j];
        int neg; iss >> neg;
        r.rat_neg = (neg != 0);
        if (iss) rels.push_back(std::move(r));
    }
    return true;
}

// ============================================================
// 並列篩 (b 方向分割)
// ============================================================

std::vector<GNFSRel> gnfsCollectParallel(
    const GNFSPoly& poly, const std::vector<uint64_t>& rfb,
    const std::vector<AlgFBEntry>& afb, const GNFSParams& par,
    int nthreads = 0)
{
    if (nthreads <= 0)
        nthreads = std::max(1, (int)std::thread::hardware_concurrency());
    if (nthreads == 1)
        return gnfsCollect(poly, rfb, afb, par);

    int rn = (int)rfb.size(), an = (int)afb.size();
    int target = rn + an + 10;

    // b 範囲をスレッドに分割
    std::vector<std::vector<GNFSRel>> thread_rels(nthreads);
    std::vector<std::vector<GNFSPartial>> thread_partials(nthreads);

    auto worker = [&](int tid) {
        uint64_t lp_bound = rfb.back() * rfb.back();
        uint64_t lp2_bound = lp_bound;
        uint64_t lp2_cofactor_bound = (lp2_bound <= 0xFFFFFFFFULL)
            ? lp2_bound * lp2_bound : 0xFFFFFFFFFFFFFFFFULL;
        int A = par.sieve_a, slen = 2 * A + 1;

        std::vector<uint64_t> mmodp(rn);
        for (int i = 0; i < rn; i++) {
            Int mm = poly.m % Int(rfb[i]);
            if (mm.isNegative()) mm += Int(rfb[i]);
            mmodp[i] = mm.toUInt64();
        }
        std::vector<float> rlog(rn), alog(an);
        for (int i = 0; i < rn; i++) rlog[i] = (float)std::log2((double)rfb[i]);
        for (int i = 0; i < an; i++) alog[i] = (float)std::log2((double)afb[i].p);

        constexpr int BLK = 8192;
        std::vector<float> sr(BLK), sa(BLK);
        auto& rels = thread_rels[tid];
        auto& parts = thread_partials[tid];

        for (int64_t b = 1 + tid; b <= par.sieve_b; b += nthreads) {
            std::vector<int64_t> rstart(rn), astart(an);
            for (int i = 0; i < rn; i++) {
                uint64_t p = rfb[i];
                uint64_t bmodp = (uint64_t)(b % (int64_t)p);
                uint64_t start = (p - bmodp * mmodp[i] % p) % p;
                rstart[i] = (int64_t)((start + (uint64_t)(A % (int64_t)p)) % p);
            }
            for (int i = 0; i < an; i++) {
                uint64_t p = afb[i].p, r = afb[i].r;
                uint64_t bmodp = (uint64_t)(b % (int64_t)p);
                uint64_t start = (p - bmodp * r % p) % p;
                astart[i] = (int64_t)((start + (uint64_t)(A % (int64_t)p)) % p);
            }

            double bm = std::abs((double)b * poly.m.toDouble());
            float rthr = (float)(std::log2(std::max(bm, 1.0) + A) - par.thresh_adj);
            float athr = (float)((double)poly.d * std::log2(std::max((double)A, 1.0)) - par.thresh_adj);
            if (rthr < 3) rthr = 3;
            if (athr < 3) athr = 3;

            for (int blk_start = 0; blk_start < slen; blk_start += BLK) {
                int blk_end = std::min(blk_start + BLK, slen);
                int blk_len = blk_end - blk_start;
                std::fill(sr.begin(), sr.begin() + blk_len, 0.0f);
                std::fill(sa.begin(), sa.begin() + blk_len, 0.0f);

                for (int i = 0; i < rn; i++) {
                    uint64_t p = rfb[i];
                    int64_t j = rstart[i];
                    if (j < blk_start) {
                        int64_t skip = (blk_start - j + (int64_t)p - 1) / (int64_t)p;
                        j += skip * (int64_t)p;
                    }
                    for (; j < blk_end; j += (int64_t)p)
                        sr[j - blk_start] += rlog[i];
                }
                for (int i = 0; i < an; i++) {
                    uint64_t p = afb[i].p;
                    int64_t j = astart[i];
                    if (j < blk_start) {
                        int64_t skip = (blk_start - j + (int64_t)p - 1) / (int64_t)p;
                        j += skip * (int64_t)p;
                    }
                    for (; j < blk_end; j += (int64_t)p)
                        sa[j - blk_start] += alog[i];
                }

                for (int il = 0; il < blk_len; il++) {
                    if (sr[il] < rthr || sa[il] < athr) continue;
                    int idx = blk_start + il;
                    int64_t a = idx - A;
                    if (a == 0) continue;
                    if (std::gcd(std::abs(a), b) != 1) continue;

                    Int rv = Int(a) + Int(b) * poly.m;
                    bool rneg = rv.isNegative();
                    Int ra = rneg ? -rv : rv;
                    std::vector<int> re(rn, 0);
                    bool r_smooth;
                    uint64_t r_rem64 = 0;
                    if (ra.bitLength() <= 64) {
                        r_rem64 = ra.toUInt64();
                        for (int i = 0; i < rn && r_rem64 > 1; i++) {
                            uint64_t p = rfb[i];
                            while (r_rem64 % p == 0) { r_rem64 /= p; re[i]++; }
                        }
                        r_smooth = (r_rem64 == 1);
                    } else {
                        Int rrem = ra;
                        for (int i = 0; i < rn; i++) {
                            Int p(rfb[i]);
                            while ((rrem % p).isZero()) { rrem /= p; re[i]++; }
                        }
                        r_smooth = (rrem == 1);
                        if (!r_smooth && rrem.bitLength() <= 64) r_rem64 = rrem.toUInt64();
                    }

                    int r_lp_count = 0;
                    uint64_t r_lps[2] = {0, 0};
                    if (!r_smooth) {
                        if (r_rem64 <= 1) continue;
                        if (r_rem64 <= lp_bound) {
                            r_lp_count = 1; r_lps[0] = r_rem64;
                        } else if (r_rem64 <= lp2_cofactor_bound) {
                            if (isPrime64(r_rem64)) continue;
                            uint64_t f1 = pollardRho64(r_rem64);
                            if (f1 == 0) continue;
                            uint64_t f2 = r_rem64 / f1;
                            if (f1 > f2) std::swap(f1, f2);
                            if (f2 > lp2_bound) continue;
                            r_lp_count = 2; r_lps[0] = f1; r_lps[1] = f2;
                        } else { continue; }
                    }

                    Int norm = algNorm(a, b, poly.f, poly.d);
                    if (norm.isZero()) continue;
                    Int na = norm.isNegative() ? -norm : norm;
                    std::vector<int> ae(an, 0);
                    bool a_smooth; uint64_t a_r64 = 0;
                    if (na.bitLength() <= 64) {
                        a_r64 = na.toUInt64();
                        for (int i = 0; i < an && a_r64 > 1; i++) {
                            uint64_t p = afb[i].p;
                            while (a_r64 % p == 0) { a_r64 /= p; ae[i]++; }
                        }
                        a_smooth = (a_r64 == 1);
                    } else {
                        Int arem = na;
                        for (int i = 0; i < an; i++) {
                            Int p(afb[i].p);
                            while ((arem % p).isZero()) { arem /= p; ae[i]++; }
                        }
                        a_smooth = (arem == 1);
                        if (!a_smooth && arem.bitLength() <= 64) a_r64 = arem.toUInt64();
                    }

                    uint64_t a_lp = 0;
                    if (!a_smooth) {
                        if (r_lp_count > 0) continue;
                        if (a_r64 <= 1 || a_r64 > lp_bound) continue;
                        a_lp = a_r64;
                    }

                    if (r_lp_count == 0 && a_lp == 0) {
                        rels.push_back({a, b, re, ae, rneg});
                    } else {
                        GNFSPartial part;
                        part.a = a; part.b = b;
                        part.rexp = re; part.aexp = ae;
                        part.rat_neg = rneg;
                        part.rat_lp[0] = r_lps[0]; part.rat_lp[1] = r_lps[1];
                        part.rat_lp_count = r_lp_count;
                        part.alg_lp = a_lp;
                        parts.push_back(std::move(part));
                    }
                }
            }
        }
    };

    // スレッド起動
    std::vector<std::thread> threads;
    for (int t = 0; t < nthreads; t++)
        threads.emplace_back(worker, t);
    for (auto& t : threads) t.join();

    // 結果をマージ
    std::vector<GNFSRel> rels;
    std::vector<GNFSPartial> partials;
    for (int t = 0; t < nthreads; t++) {
        rels.insert(rels.end(), thread_rels[t].begin(), thread_rels[t].end());
        partials.insert(partials.end(), thread_partials[t].begin(), thread_partials[t].end());
    }

    // LP マージ (gnfsCollect と同じ 3 フェーズ)
    // Phase 1: Rational 1-LP
    {
        std::unordered_map<uint64_t, size_t> lp_map;
        for (size_t i = 0; i < partials.size(); i++) {
            if (partials[i].rat_lp_count != 1 || partials[i].alg_lp != 0) continue;
            auto it = lp_map.find(partials[i].rat_lp[0]);
            if (it != lp_map.end()) {
                auto& p1 = partials[it->second]; auto& p2 = partials[i];
                GNFSRel c; c.a = p1.a; c.b = p1.b;
                c.rexp.resize(rn, 0); c.aexp.resize(an, 0);
                for (int j = 0; j < rn; j++) c.rexp[j] = p1.rexp[j] + p2.rexp[j];
                for (int j = 0; j < an; j++) c.aexp[j] = p1.aexp[j] + p2.aexp[j];
                c.rat_neg = (p1.rat_neg != p2.rat_neg);
                rels.push_back(c);
            } else lp_map[partials[i].rat_lp[0]] = i;
        }
    }
    // Phase 2: Algebraic 1-LP
    {
        std::unordered_map<uint64_t, size_t> alp_map;
        for (size_t i = 0; i < partials.size(); i++) {
            if (partials[i].rat_lp_count != 0 || partials[i].alg_lp == 0) continue;
            auto it = alp_map.find(partials[i].alg_lp);
            if (it != alp_map.end()) {
                auto& p1 = partials[it->second]; auto& p2 = partials[i];
                GNFSRel c; c.a = p1.a; c.b = p1.b;
                c.rexp.resize(rn, 0); c.aexp.resize(an, 0);
                for (int j = 0; j < rn; j++) c.rexp[j] = p1.rexp[j] + p2.rexp[j];
                for (int j = 0; j < an; j++) c.aexp[j] = p1.aexp[j] + p2.aexp[j];
                c.rat_neg = (p1.rat_neg != p2.rat_neg);
                rels.push_back(c);
            } else alp_map[partials[i].alg_lp] = i;
        }
    }

    return rels;
}

// ============================================================
// Algebraic square root via Hensel lifting
// ============================================================

// splitting prime を探す (f が d 個の異なる根を持つ p)
// f の判別式 disc(f) ≢ 0 (mod p) かつ f が d 個の異なる根を持つ素数 p を返す
uint64_t findSplitPrime(const Poly& f, int d) {
    // エラトステネスの篩で候補素数を生成
    constexpr uint64_t SEARCH_LIMIT = 100000;
    std::vector<bool> is_prime(SEARCH_LIMIT + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (uint64_t i = 2; i * i <= SEARCH_LIMIT; i++)
        if (is_prime[i])
            for (uint64_t j = i * i; j <= SEARCH_LIMIT; j += i)
                is_prime[j] = false;

    // 主係数 f[d] mod p ≠ 0 かつ完全分裂する素数を探す
    for (uint64_t p = 101; p <= SEARCH_LIMIT; p++) {
        if (!is_prime[p]) continue;
        // f[d] mod p ≠ 0 を確認 (判別式の必要条件)
        Int fd_mod = f[d] % Int(p);
        if (fd_mod.isZero()) continue;
        auto roots = findRootsModP(f, d, p);
        if ((int)roots.size() == d) return p;
    }
    return 0;
}

std::vector<Int> algSqrtCandidates(
    const std::vector<std::pair<int64_t,int64_t>>& pairs,
    const Poly& f, int d, const Int& m, const Int& n)
{
    uint64_t p = findSplitPrime(f, d);
    if (p == 0) return {};

    // f の根 mod p (Cantor-Zassenhaus)
    auto roots = findRootsModP(f, d, p);
    if ((int)roots.size() != d) return {};

    // 目標 modulus: p^k > 推定係数上限
    // 簡易上限: n^{pairs.size()/d + 1}
    int R = (int)pairs.size();
    Int pk = Int(p);
    // p^k > n * 10^(R) くらいを目標 (保守的)
    Int target_pk = n * n; // 最低でも n^2
    // 実用的上限: R * log2(max(a,b)) bits
    int lift_bits = std::max(200, R * 15 / d);
    while ((int)pk.bitLength() < lift_bits) pk = pk * pk;

    // Π(x) = ∏(a_i + b_i·x) mod f(x) mod pk を計算
    Poly Pi(d, Int(0));
    Pi[0] = 1;
    for (auto& [a, b] : pairs) {
        Poly fac(d, Int(0));
        fac[0] = Int(a); if (fac[0].isNegative()) fac[0] = ((fac[0] % pk) + pk) % pk;
        if (d > 1) {
            fac[1] = Int(b);
            if (fac[1].isNegative()) fac[1] = ((fac[1] % pk) + pk) % pk;
        }
        Pi = polyMul(Pi, fac, f, d, pk);
    }

    // Π(root_j) mod p の計算
    Int pmod(p);
    std::vector<uint64_t> pi_at_r(d);
    for (int j = 0; j < d; j++) {
        uint64_t val = 0, rp = 1;
        for (int i = 0; i < d; i++) {
            val = (val + (Pi[i] % pmod).toUInt64() * rp) % p;
            rp = rp * roots[j] % p;
        }
        pi_at_r[j] = val;
    }

    // sqrt(Π(root_j)) mod p
    std::vector<uint64_t> sqrts(d);
    for (int j = 0; j < d; j++) {
        sqrts[j] = gnfs_sqrtmod(pi_at_r[j], p);
        if (sqrts[j] * sqrts[j] % p != pi_at_r[j]) return {}; // 非 QR
    }

    // 2^d 通りの符号組み合わせを試す
    std::vector<Int> candidates;
    int ncomb = 1 << d;

    for (int combo = 0; combo < ncomb; combo++) {
        // Lagrange 補間で h(x) mod p を構成
        std::vector<uint64_t> ys(d);
        for (int j = 0; j < d; j++)
            ys[j] = (combo & (1 << j)) ? (p - sqrts[j]) % p : sqrts[j];

        Poly h(d, Int(0));
        for (int j = 0; j < d; j++) {
            uint64_t denom = 1;
            for (int k = 0; k < d; k++) {
                if (k == j) continue;
                denom = denom * ((roots[j] + p - roots[k]) % p) % p;
            }
            uint64_t scale = ys[j] * gnfs_modinv(denom, p) % p;

            // L_j(x) = ∏_{k≠j}(x - roots[k]) / denom
            Poly Lj(d, Int(0)); Lj[0] = 1; int deg = 0;
            for (int k = 0; k < d; k++) {
                if (k == j) continue;
                Poly nLj(d, Int(0));
                for (int c = deg+1; c >= 1; c--) nLj[c] = Lj[c-1];
                uint64_t neg_rk = (p - roots[k]) % p;
                for (int c = 0; c <= deg; c++)
                    nLj[c] = (nLj[c] + Lj[c] * Int(neg_rk)) % pmod;
                Lj = nLj; deg++;
            }
            for (int c = 0; c < d; c++)
                h[c] = (h[c] + Lj[c] * Int(scale)) % pmod;
        }
        for (int c = 0; c < d; c++) if (h[c].isNegative()) h[c] += pmod;

        // (2h)^{-1} mod f mod p を補間
        Poly u(d, Int(0));
        for (int j = 0; j < d; j++) {
            uint64_t hval = 0, rp = 1;
            for (int i = 0; i < d; i++) {
                hval = (hval + h[i].toUInt64() * rp) % p;
                rp = rp * roots[j] % p;
            }
            uint64_t inv2h = gnfs_modinv(2 * hval % p, p);

            uint64_t denom = 1;
            for (int k = 0; k < d; k++) {
                if (k == j) continue;
                denom = denom * ((roots[j] + p - roots[k]) % p) % p;
            }
            uint64_t scale = inv2h * gnfs_modinv(denom, p) % p;

            Poly Lj(d, Int(0)); Lj[0] = 1; int deg = 0;
            for (int k = 0; k < d; k++) {
                if (k == j) continue;
                Poly nLj(d, Int(0));
                for (int c = deg+1; c >= 1; c--) nLj[c] = Lj[c-1];
                uint64_t neg_rk = (p - roots[k]) % p;
                for (int c = 0; c <= deg; c++)
                    nLj[c] = (nLj[c] + Lj[c] * Int(neg_rk)) % pmod;
                Lj = nLj; deg++;
            }
            for (int c = 0; c < d; c++)
                u[c] = (u[c] + Lj[c] * Int(scale)) % pmod;
        }
        for (int c = 0; c < d; c++) if (u[c].isNegative()) u[c] += pmod;

        // Hensel lifting: h, u を mod p から mod pk へ
        Int cur_mod = pmod;
        while (cur_mod < pk) {
            Int new_mod = cur_mod * cur_mod;
            if (new_mod > pk) new_mod = pk;

            // Π mod new_mod
            Poly Pi_nm(d);
            for (int i = 0; i < d; i++) Pi_nm[i] = Pi[i] % new_mod;

            // h' = h - (h² - Π)·u mod f mod new_mod
            Poly h2 = polyMul(h, h, f, d, new_mod);
            Poly diff(d);
            for (int i = 0; i < d; i++) {
                diff[i] = (h2[i] - Pi_nm[i]) % new_mod;
                if (diff[i].isNegative()) diff[i] += new_mod;
            }
            Poly corr = polyMul(diff, u, f, d, new_mod);
            for (int i = 0; i < d; i++) {
                h[i] = (h[i] - corr[i]) % new_mod;
                if (h[i].isNegative()) h[i] += new_mod;
            }

            // u' = u·(2 - 2h·u) mod f mod new_mod
            Poly two_h(d);
            for (int i = 0; i < d; i++) two_h[i] = (2 * h[i]) % new_mod;
            Poly hu = polyMul(two_h, u, f, d, new_mod);
            Poly tmhu(d);
            tmhu[0] = (2 - hu[0]) % new_mod;
            if (tmhu[0].isNegative()) tmhu[0] += new_mod;
            for (int i = 1; i < d; i++) {
                tmhu[i] = new_mod - hu[i];
                if (tmhu[i] == new_mod) tmhu[i] = Int(0);
            }
            u = polyMul(u, tmhu, f, d, new_mod);

            cur_mod = new_mod;
        }

        // Balanced representation [-pk/2, pk/2]
        Int half_pk = pk / 2;
        for (int i = 0; i < d; i++)
            if (h[i] > half_pk) h[i] -= pk;

        // Y = h(m) mod n
        Int Y(0), mpow(1);
        for (int i = 0; i < d; i++) {
            Y = (Y + h[i] * mpow) % n;
            mpow = (mpow * m) % n;
        }
        if (Y.isNegative()) Y += n;
        candidates.push_back(Y);
    }
    return candidates;
}

} // anonymous namespace

// ============================================================
// findFactor_GNFS — GNFS メインエントリポイント
// ============================================================

Int IntPrime::findFactor_GNFS(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;

    std::string ns = n.toString();
    int digits = (int)ns.size();
    if (digits < 25) return findFactor_SIQS(n);

    auto par = gnfsParams(digits);

    // 1. Polynomial selection (Kleinjung lattice-based for 40+, baseMSelect for small)
    auto poly = (digits >= 40) ? kleinjungSelect(n, par.degree)
                               : baseMSelect(n, par.degree);

    // Verify f(m) = n
    {
        Int check(0), mp(1);
        for (int i = 0; i <= poly.d; i++) { check += poly.f[i] * mp; mp *= poly.m; }
        if (check != n) return findFactor_SIQS(n);
    }

    // 2. Factor bases
    auto rfb = buildRatFB(par.rfb_bound);
    auto afb = buildAlgFB(poly.f, poly.d, rfb);
    int rn = (int)rfb.size(), an = (int)afb.size();
    int total_cols = rn + an + 2;

    // 3. Relation collection (並列篩 → 不足なら lattice sieve で補完)
    auto rels = gnfsCollectParallel(poly, rfb, afb, par);
    if ((int)rels.size() < total_cols + 1) {
        // Lattice sieve で追加関係を収集
        std::vector<GNFSPartial> extra_partials;
        uint64_t lp_bound = rfb.back() * rfb.back();
        gnfsLatticeSieve(poly, rfb, afb, par, rels, extra_partials,
                         lp_bound, total_cols + 10);
        // 追加 partial の 1-LP ペアリング
        {
            std::unordered_map<uint64_t, size_t> lp_map;
            for (size_t i = 0; i < extra_partials.size(); i++) {
                if (extra_partials[i].rat_lp_count != 1) continue;
                auto it = lp_map.find(extra_partials[i].rat_lp[0]);
                if (it != lp_map.end()) {
                    auto& p1 = extra_partials[it->second];
                    auto& p2 = extra_partials[i];
                    GNFSRel c; c.a = p1.a; c.b = p1.b;
                    c.rexp.resize(rn, 0); c.aexp.resize(an, 0);
                    for (int j = 0; j < rn; j++) c.rexp[j] = p1.rexp[j] + p2.rexp[j];
                    for (int j = 0; j < an; j++) c.aexp[j] = p1.aexp[j] + p2.aexp[j];
                    c.rat_neg = (p1.rat_neg != p2.rat_neg);
                    rels.push_back(c);
                } else lp_map[extra_partials[i].rat_lp[0]] = i;
            }
        }
    }
    if ((int)rels.size() < total_cols + 1)
        return findFactor_SIQS(n);

    // 4. Linear algebra GF(2)
    int nrels = (int)rels.size();
    std::vector<std::vector<uint8_t>> mat(nrels, std::vector<uint8_t>(total_cols, 0));
    for (int i = 0; i < nrels; i++) {
        for (int j = 0; j < rn; j++) mat[i][j] = rels[i].rexp[j] & 1;
        for (int j = 0; j < an; j++) mat[i][rn + j] = rels[i].aexp[j] & 1;
        mat[i][rn + an] = rels[i].rat_neg ? 1 : 0;
        // algebraic sign column
        Int norm = algNorm(rels[i].a, rels[i].b, poly.f, poly.d);
        mat[i][rn + an + 1] = norm.isNegative() ? 1 : 0;
    }

    auto null_vecs = gf2_linalg::findNullSpaceGF2(mat, nrels, total_cols);

    // 5. Factor extraction
    for (auto& indices : null_vecs) {
        // Rational square root: X = ∏ p^{e/2} mod n
        std::vector<int> tr(rn, 0);
        for (int idx : indices)
            for (int j = 0; j < rn; j++) tr[j] += rels[idx].rexp[j];

        Int X(1);
        for (int j = 0; j < rn; j++) {
            int half = tr[j] / 2;
            if (half > 0) X = (X * IntModular::powerMod(Int(rfb[j]), Int(half), n)) % n;
        }

        // Algebraic square root
        std::vector<std::pair<int64_t,int64_t>> ab;
        for (int idx : indices) ab.push_back({rels[idx].a, rels[idx].b});
        auto Ys = algSqrtCandidates(ab, poly.f, poly.d, poly.m, n);

        for (auto& Y : Ys) {
            Int g = IntGCD::gcd((X - Y + n) % n, n);
            if (g > 1 && g < n) return g;
            g = IntGCD::gcd((X + Y) % n, n);
            if (g > 1 && g < n) return g;
        }
    }

    return findFactor_SIQS(n); // Fallback
}

// ============================================================
// 分散篩インフラ (Distributed Sieving Infrastructure)
// ============================================================
// 篩タスクを b 範囲で分割し、独立プロセスで実行可能にする。
// 各ワーカーは saveSieveTask → sieve → saveSieveResult のサイクル。
// マスターは mergeSieveResults で統合。

namespace distributed {

/// 篩タスクの仕様 (ワーカーへ配布)
struct SieveTask {
    Int n;                  // 因数分解対象
    GNFSPoly poly;          // 多項式
    int b_start, b_end;     // 担当 b 範囲 [b_start, b_end)
    int sieve_a;            // a 範囲 [-sieve_a, sieve_a]
    int rfb_bound;          // 有理 FB 上限
    int thresh_adj;         // 篩閾値調整
};

/// 篩タスクをファイルに保存
bool saveSieveTask(const std::string& path, const SieveTask& task) {
    std::ofstream ofs(path);
    if (!ofs) return false;
    ofs << "GNFS_SIEVE_TASK_V1\n";
    ofs << "n " << task.n.toString() << "\n";
    ofs << "degree " << task.poly.d << "\n";
    ofs << "poly";
    for (int i = 0; i <= task.poly.d; i++)
        ofs << " " << task.poly.f[i].toString();
    ofs << "\n";
    ofs << "m " << task.poly.m.toString() << "\n";
    ofs << "b_range " << task.b_start << " " << task.b_end << "\n";
    ofs << "sieve_a " << task.sieve_a << "\n";
    ofs << "rfb_bound " << task.rfb_bound << "\n";
    ofs << "thresh_adj " << task.thresh_adj << "\n";
    return ofs.good();
}

/// 篩タスクをファイルから読み込み
bool loadSieveTask(const std::string& path, SieveTask& task) {
    std::ifstream ifs(path);
    if (!ifs) return false;
    std::string line;
    if (!std::getline(ifs, line) || line != "GNFS_SIEVE_TASK_V1")
        return false;

    while (std::getline(ifs, line)) {
        std::istringstream iss(line);
        std::string key;
        iss >> key;
        if (key == "n") {
            std::string s; iss >> s;
            task.n = Int(s);
        } else if (key == "degree") {
            iss >> task.poly.d;
        } else if (key == "poly") {
            task.poly.f.resize(task.poly.d + 1);
            for (int i = 0; i <= task.poly.d; i++) {
                std::string s; iss >> s;
                task.poly.f[i] = Int(s);
            }
        } else if (key == "m") {
            std::string s; iss >> s;
            task.poly.m = Int(s);
        } else if (key == "b_range") {
            iss >> task.b_start >> task.b_end;
        } else if (key == "sieve_a") {
            iss >> task.sieve_a;
        } else if (key == "rfb_bound") {
            iss >> task.rfb_bound;
        } else if (key == "thresh_adj") {
            iss >> task.thresh_adj;
        }
    }
    return true;
}

/// 篩の b 範囲を num_workers 分割してタスクリストを生成
std::vector<SieveTask> splitSieveTasks(
    const Int& n, const GNFSPoly& poly, const GNFSParams& par,
    int num_workers)
{
    std::vector<SieveTask> tasks;
    int total_b = par.sieve_b;
    int chunk = std::max(1, total_b / num_workers);

    for (int w = 0; w < num_workers; w++) {
        SieveTask t;
        t.n = n;
        t.poly = poly;
        t.b_start = w * chunk + 1;
        t.b_end = (w == num_workers - 1) ? total_b + 1 : (w + 1) * chunk + 1;
        t.sieve_a = par.sieve_a;
        t.rfb_bound = par.rfb_bound;
        t.thresh_adj = par.thresh_adj;
        tasks.push_back(t);
    }
    return tasks;
}

/// 複数の関係ファイルを統合
bool mergeSieveResults(
    const std::vector<std::string>& paths,
    int rn, int an,
    std::vector<GNFSRel>& merged_rels)
{
    for (auto& p : paths) {
        if (!gnfsLoadRelations(p, rn, an, merged_rels))
            return false;
    }
    return true;
}

} // namespace distributed

} // namespace calx