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

// MultivariatePolynomial_factorization.hpp
// 多変数多項式の因数分解 (Wang EEZ + Hensel lifting)
//
// CAS-3b-1: Hensel lifting (多変数版)
// CAS-3b-2: 評価点選択 (Leading Coefficient Problem)
//
// アルゴリズム:
//   1. content-primitive 分離
//   2. 評価点選択 (補助変数に小整数を代入して 1 変数化)
//   3. 1 変数因数分解 (既存の factorize を利用)
//   4. Hensel lifting で多変数因子を復元
//
// 参考文献:
//   Wang "An Improved Multivariate Polynomial Factoring Algorithm" (1978)
//   Geddes, Czapor, Labahn "Algorithms for Computer Algebra" Ch.6

#pragma once

#include <math/core/MultivariatePolynomial.hpp>
#include <math/core/Polynomial.hpp>
#include <math/core/Matrix.hpp>
#include <math/core/Polynomial_factorization.hpp>
#include <math/core/FiniteField.hpp>

#include <cstdint>
#include <vector>
#include <set>
#include <algorithm>
#include <numeric>
#include <functional>

namespace calx {

namespace detail {

// Monomial 構築ヘルパー
inline Monomial makeSingleVarMono(int var, int exp) {
    std::vector<int> e(var + 1, 0);
    e[var] = exp;
    return Monomial(std::move(e));
}

// 多変数多項式を 1 変数多項式に変換 (主変数のみ残し、補助変数に値を代入)
template<typename T>
Polynomial<int64_t> mpolyToUnivariate(
    const MultivariatePolynomial<T>& f, int mainVar,
    const std::vector<std::pair<int, T>>& evalPoints)
{
    auto reduced = f;
    for (auto& [var, val] : evalPoints)
        reduced = reduced.substitute(var, val);

    // mainVar の 1 変数多項式に変換
    auto coeffs = viewAsUnivariate(reduced, mainVar);
    std::vector<int64_t> uni_c;
    for (auto& c : coeffs) {
        if (c.isZero()) uni_c.push_back(0);
        else uni_c.push_back(static_cast<int64_t>(c.terms()[0].coeff));
    }
    while (uni_c.size() > 1 && uni_c.back() == 0) uni_c.pop_back();
    return Polynomial<int64_t>(uni_c);
}

// 1 変数因子を多変数多項式に変換
template<typename T>
MultivariatePolynomial<T> uniToMpoly(const Polynomial<int64_t>& f, int mainVar)
{
    MultivariatePolynomial<T> result;
    for (int i = 0; i <= f.degree(); ++i) {
        if (f[i] == 0) continue;
        result.addTerm(makeSingleVarMono(mainVar, i), static_cast<T>(f[i]));
    }
    return result;
}

// 二変数 Hensel lifting
// f(x,y) ≡ g₀(x)·h₀(x) mod (y - a) を mod (y - a)^{d+1} まで持ち上げ
// g₀, h₀ は y=a での 1 変数因子
// 出力: g(x,y), h(x,y) s.t. f = g·h
// 二変数 Hensel lifting
// 前提: f(x, y) は y=0 で g0·h0 に分解済み (呼び出し側で変数シフト済み)
template<typename T>
bool bivarHenselLift(
    const MultivariatePolynomial<T>& f,
    MultivariatePolynomial<T>& g_out, MultivariatePolynomial<T>& h_out,
    const Polynomial<int64_t>& g0, const Polynomial<int64_t>& h0,
    int mainVar, int auxVar)
{
    using MPoly = MultivariatePolynomial<T>;
    // f の auxVar に関する最大次数
    int maxDeg = 0;
    for (const auto& t : f.terms())
        maxDeg = std::max(maxDeg, t.mono[auxVar]);

    // f を auxVar の冪で展開: f = Σ f_j(x) · y^j
    // (呼び出し側で変数シフト済みなので auxVal=0 前提)
    auto f_coeffs = viewAsUnivariate(f, auxVar);

    // g, h を y の冪で表現
    // g = g₀(x) + g₁(x)·y + g₂(x)·y² + ...
    // h = h₀(x) + h₁(x)·y + h₂(x)·y² + ...
    // 初期値 (y=a で): g₀ と h₀ は既に持っている

    // 1 変数の係数列
    int gd = g0.degree(), hd = h0.degree();
    // g_coeffs[j] = x の多項式 (g の y^j 次の係数)
    std::vector<Polynomial<int64_t>> gc(maxDeg + 1, Polynomial<int64_t>({0}));
    std::vector<Polynomial<int64_t>> hc(maxDeg + 1, Polynomial<int64_t>({0}));
    gc[0] = g0;
    hc[0] = h0;

    // 注意: gc, hc の全係数は mod p_work で管理する
    // 畳み込み: (g·h)_j = Σ_{k=0}^{j} g_k · h_{j-k}
    // f_j = (g·h)_j が成立するように g_j, h_j を決定

    // 素数の選択: 結果の係数が対称剰余 [-p/2, p/2] 内に収まるよう
    // f の最大係数を元に bound を推定
    int64_t max_coeff = 1;
    for (const auto& t : f.terms()) {
        int64_t c = static_cast<int64_t>(t.coeff);
        if (c < 0) c = -c;
        max_coeff = std::max(max_coeff, c);
    }
    int64_t bound = max_coeff * int64_t(maxDeg + 1) * int64_t(g0.degree() + h0.degree() + 2);
    int64_t p_work = std::max(bound * 4 + 1, int64_t(1009));
    // 素数に繰り上げ
    auto isPrimeSmall = [](int64_t n) -> bool {
        if (n < 2) return false;
        if (n < 4) return true;
        if (n % 2 == 0 || n % 3 == 0) return false;
        for (int64_t i = 5; i * i <= n; i += 6)
            if (n % i == 0 || n % (i + 2) == 0) return false;
        return true;
    };
    while (!isPrimeSmall(p_work)) ++p_work;
    auto g0_mod = g0, h0_mod = h0;

    // Extended GCD for Bezout
    Polynomial<int64_t> gcd_p, s_p, t_p;

    // Extended GCD mod p_work: s·g₀ + t·h₀ ≡ 1 (mod p_work)
    Polynomial<int64_t> s_mod, t_mod;
    {
        auto reduceMod = [p_work](Polynomial<int64_t>& f) {
            std::vector<int64_t> c;
            for (int i = 0; i <= f.degree(); ++i) {
                int64_t v = f[i] % p_work;
                if (v < 0) v += p_work;
                if (v > p_work / 2) v -= p_work;
                c.push_back(v);
            }
            while (c.size() > 1 && c.back() == 0) c.pop_back();
            f = Polynomial<int64_t>(c);
        };
        auto mulMod = [p_work](const Polynomial<int64_t>& a, const Polynomial<int64_t>& b) {
            auto r = a * b;
            std::vector<int64_t> c;
            for (int i = 0; i <= r.degree(); ++i) {
                int64_t v = r[i] % p_work;
                if (v < 0) v += p_work;
                if (v > p_work / 2) v -= p_work;
                c.push_back(v);
            }
            while (c.size() > 1 && c.back() == 0) c.pop_back();
            return Polynomial<int64_t>(c);
        };
        auto modInv = [p_work](int64_t a) -> int64_t {
            a = ((a % p_work) + p_work) % p_work;
            int64_t r0 = a, r1 = p_work, s0 = 1, s1 = 0;
            while (r1 != 0) {
                int64_t q = r0 / r1, tmp;
                tmp = r1; r1 = r0 - q * r1; r0 = tmp;
                tmp = s1; s1 = s0 - q * s1; s0 = tmp;
            }
            return ((s0 % p_work) + p_work) % p_work;
        };
        auto divMod = [&](const Polynomial<int64_t>& a, const Polynomial<int64_t>& b)
            -> std::pair<Polynomial<int64_t>, Polynomial<int64_t>> {
            if (a.degree() < b.degree()) return {Polynomial<int64_t>({0}), a};
            std::vector<int64_t> rem(a.degree() + 1);
            for (int i = 0; i <= a.degree(); ++i) rem[i] = ((a[i] % p_work) + p_work) % p_work;
            int64_t lc_inv = modInv(b[b.degree()]);
            int dq = a.degree() - b.degree();
            std::vector<int64_t> q(dq + 1, 0);
            for (int i = dq; i >= 0; --i) {
                q[i] = (rem[i + b.degree()] * lc_inv) % p_work;
                for (int j = 0; j <= b.degree(); ++j)
                    rem[i + j] = ((rem[i + j] - q[i] * (b[j] % p_work)) % p_work + p_work) % p_work;
            }
            int64_t half = p_work / 2;
            for (auto& v : q) { if (v > half) v -= p_work; }
            for (auto& v : rem) { if (v > half) v -= p_work; }
            while (q.size() > 1 && q.back() == 0) q.pop_back();
            while (rem.size() > 1 && rem.back() == 0) rem.pop_back();
            return {Polynomial<int64_t>(q), Polynomial<int64_t>(rem)};
        };

        Polynomial<int64_t> r0 = g0, r1 = h0;
        Polynomial<int64_t> s0({1}), s1({0});
        Polynomial<int64_t> t0({0}), t1({1});
        reduceMod(r0); reduceMod(r1);
        for (int iter = 0; iter < 100; ++iter) {
            bool rz = true;
            for (int i = 0; i <= r1.degree(); ++i) if (r1[i] != 0) { rz = false; break; }
            if (rz) break;
            auto [q, r] = divMod(r0, r1);
            r0 = r1; r1 = r;
            auto ns = s0 - mulMod(q, s1); reduceMod(ns); s0 = s1; s1 = ns;
            auto nt = t0 - mulMod(q, t1); reduceMod(nt); t0 = t1; t1 = nt;
        }
        // gcd を 1 に正規化: s·g0 + t·h0 = gcd → (s/gcd)·g0 + (t/gcd)·h0 = 1
        // gcd は定数多項式のはず (g0, h0 が coprime)
        if (!r0.isZero() && r0[0] != 0) {
            int64_t gcd_inv = modInv(((r0[0] % p_work) + p_work) % p_work);
            auto scalePoly = [&](Polynomial<int64_t>& f) {
                std::vector<int64_t> c;
                for (int i = 0; i <= f.degree(); ++i) {
                    int64_t v = (f[i] * gcd_inv) % p_work;
                    if (v < 0) v += p_work;
                    if (v > p_work / 2) v -= p_work;
                    c.push_back(v);
                }
                while (c.size() > 1 && c.back() == 0) c.pop_back();
                f = Polynomial<int64_t>(c);
            };
            scalePoly(s0);
            scalePoly(t0);
        }
        s_mod = s0;
        t_mod = t0;

        // Bezout 分配: g0·h_j + h0·g_j = delta
        // g_j = (s·delta) mod h0 (mod p_work)
        // h_j = (t·delta) mod g0 (mod p_work)
        // 結果を対称剰余に戻してから、真の小さい整数係数を復元
        // distribute: δ = g0·h_j + h0·g_j を解く
        // s·g0 + t·h0 = 1 (mod p) より:
        //   g_j = (s·δ) rem h0  (deg < deg(h0))
        //   h_j = (t·δ + q·g0)  where q = (s·δ) quo h0
        // これにより g0·h_j + h0·g_j = δ が成立
        // s·δ = h_j + h0·(...) → (s·δ) rem h0 = h_j (deg < deg(h0))
        // t·δ = g_j + g0·(...) → (t·δ) rem g0 = g_j (deg < deg(g0))
        // 残り: g_j の補正 = t·δ + q_s·g0 (q_s = (s·δ) quo h0)
        // → 実は g_j = (t·δ + q_s·g0) が deg < deg(g0) とは限らない
        // 正しい分配:
        //   h_j = (s·δ) rem h0
        //   g_j = (t·δ + q·g0)  where q = (s·δ) quo h0
        auto distribute = [&](const Polynomial<int64_t>& delta) {
            auto sd = mulMod(s_mod, delta);
            auto [q_sd, r_sd] = divMod(sd, h0);   // h_j = r_sd
            auto td = mulMod(t_mod, delta);
            auto g_j = td + mulMod(q_sd, g0);     // g_j = t·δ + q·g0
            reduceMod(g_j);
            return std::make_pair(g_j, r_sd);     // (g_j, h_j)
        };

        for (int j = 1; j <= maxDeg; ++j) {
            // 畳み込みの既知部分 (mod p_work)
            Polynomial<int64_t> conv({0});
            for (int k = 1; k < j; ++k) {
                auto term = mulMod(gc[k], hc[j - k]);
                conv = conv + term;
                reduceMod(conv);
            }

            // f_j(x)
            Polynomial<int64_t> fj({0});
            if (j < static_cast<int>(f_coeffs.size()) && !f_coeffs[j].isZero()) {
                auto fj_u = viewAsUnivariate(f_coeffs[j], mainVar);
                std::vector<int64_t> fj_c;
                for (auto& c : fj_u) {
                    if (c.isZero()) fj_c.push_back(0);
                    else fj_c.push_back(static_cast<int64_t>(c.terms()[0].coeff));
                }
                while (fj_c.size() > 1 && fj_c.back() == 0) fj_c.pop_back();
                fj = Polynomial<int64_t>(fj_c);
            }
            reduceMod(fj);

            Polynomial<int64_t> delta = fj - conv;
            reduceMod(delta);
            auto [gj, hj] = distribute(delta);
            gc[j] = gj;
            hc[j] = hj;
        }
    }

    // g, h を多変数多項式に再構築
    // gc[j], hc[j] は (y-a)^j の係数 → y の冪に展開
    // (y-a)^j = Σ_{m=0}^{j} C(j,m)·(-a)^{j-m}·y^m
    // gc[j], hc[j] は y^j の係数 (mod p_work, auxVal=0 前提)
    // 多変数多項式に再構築 (対称剰余で整数に戻す)
    int64_t half = p_work / 2;
    auto toMPoly = [&](const std::vector<Polynomial<int64_t>>& coeffs) -> MPoly {
        MPoly result;
        for (int j = 0; j <= maxDeg; ++j) {
            if (coeffs[j].isZero()) continue;
            for (int i = 0; i <= coeffs[j].degree(); ++i) {
                int64_t c = coeffs[j][i] % p_work;
                if (c < 0) c += p_work;
                if (c > half) c -= p_work;
                if (c == 0) continue;
                std::vector<int> ev(std::max(mainVar, auxVar) + 1, 0);
                ev[mainVar] = i;
                ev[auxVar] = j;
                result.addTerm(Monomial(std::move(ev)), static_cast<T>(c));
            }
        }
        return result;
    };

    MPoly g_result = toMPoly(gc);
    MPoly h_result = toMPoly(hc);

    // 検証
    auto product = g_result * h_result;
    if (f == product) {
        g_out = g_result;
        h_out = h_result;
        return true;
    }

    // 係数補正: mod p_work で計算したため、真の係数とずれている可能性
    // 小さい係数なら正しいはず → 検証失敗は Bezout 分配の精度問題
    return false;
}

} // namespace detail

// ================================================================
// factorMultivariate メインエントリ
// ================================================================

/// 多変数多項式の因数分解 (Wang EEZ + Hensel lifting)
template<typename T>
std::vector<MultivariatePolynomial<T>> factorMultivariate(
    const MultivariatePolynomial<T>& f)
{
    static_assert(std::is_integral_v<T>,
                  "factorMultivariate requires integral coefficient type");

    using MPoly = MultivariatePolynomial<T>;

    if (f.isZero()) return {};
    if (f.terms().size() <= 1) {
        // 単項式: 各変数の冪を分離
        const auto& t = f.terms()[0];
        std::vector<MPoly> factors;
        T c = t.coeff;
        if (c != T(1)) {
            factors.push_back(MPoly(c));
        }
        auto ev = t.mono.exponents();
        for (size_t i = 0; i < ev.size(); ++i) {
            for (int e = 0; e < ev[i]; ++e) {
                factors.push_back(MPoly::variable(static_cast<int>(i)));
            }
        }
        if (factors.empty()) return {f};
        return factors;
    }

    // 使用している変数を特定
    std::set<int> vars;
    for (const auto& t : f.terms())
        for (int i = 0; i < static_cast<int>(t.mono.exponents().size()); ++i)
            if (t.mono[i] > 0) vars.insert(i);

    if (vars.empty()) return {f};

    // 整数 content を抽出
    T cont = T(0);
    for (const auto& t : f.terms()) {
        T ac = (t.coeff < T(0)) ? -t.coeff : t.coeff;
        if (cont == T(0)) cont = ac;
        else cont = std::gcd(cont, ac);
    }
    MPoly prim = f;
    if (cont > T(1)) prim = prim / cont;
    if (!prim.isZero() && prim.leadingCoefficient() < T(0)) {
        prim = prim * T(-1);
        cont = T(0) - cont;
    }

    // 1 変数のみ → 1 変数因数分解に委譲
    if (vars.size() == 1) {
        int var = *vars.begin();
        auto coeffs = viewAsUnivariate(prim, var);
        std::vector<int64_t> uniCoeffs;
        for (auto& c : coeffs) {
            if (c.isZero()) uniCoeffs.push_back(0);
            else uniCoeffs.push_back(static_cast<int64_t>(c.terms()[0].coeff));
        }
        Polynomial<int64_t> uniPoly(uniCoeffs);
        auto factored = factorize(uniPoly);

        std::vector<MPoly> factors;
        if (cont != T(1) && cont != T(0))
            factors.push_back(MPoly(cont));

        for (size_t i = 0; i < factored.size(); ++i) {
            auto& fi = factored[i];
            if (fi.degree() <= 0) {
                if (fi[0] != 0 && fi[0] != 1)
                    factors.push_back(MPoly(static_cast<T>(fi[0])));
                continue;
            }
            MPoly mf;
            for (int k = 0; k <= fi.degree(); ++k) {
                if (fi[k] == 0) continue;
                mf.addTerm(detail::makeSingleVarMono(var, k), static_cast<T>(fi[k]));
            }
            int exp_i = factored.exponent(i);
            for (int e = 0; e < exp_i; ++e)
                factors.push_back(mf);
        }
        if (factors.empty()) return {f};
        return factors;
    }

    // 多変数: 主変数選択 (最高次数)
    int mainVar = -1;
    int maxDeg = -1;
    for (int v : vars) {
        int deg = 0;
        for (const auto& t : prim.terms())
            deg = std::max(deg, t.mono[v]);
        if (deg > maxDeg) { maxDeg = deg; mainVar = v; }
    }

    // 補助変数
    std::vector<int> auxVars;
    for (int v : vars)
        if (v != mainVar) auxVars.push_back(v);

    // 評価点選択: 小整数を試す
    std::vector<T> candidates = {T(0), T(1), T(-1), T(2), T(-2), T(3), T(-3)};

    // 全補助変数に対する評価点の組み合わせを試す
    std::vector<std::pair<int, T>> bestPoints;
    Polynomial<int64_t> bestUni;
    std::vector<Polynomial<int64_t>> bestFactors;
    bool found = false;

    std::function<void(int, std::vector<std::pair<int, T>>&)>
    tryPoints = [&](int idx, std::vector<std::pair<int, T>>& points) {
        if (found) return;
        if (idx == static_cast<int>(auxVars.size())) {
            auto f_uni = detail::mpolyToUnivariate<T>(prim, mainVar, points);
            if (f_uni.degree() < maxDeg) return;  // lc が消えた

            // 無平方チェック
            if (f_uni.degree() >= 2) {
                auto df = f_uni.derivative();
                auto g = gcd(f_uni, df);
                if (g.degree() > 0) return;
            }

            auto factored = factorize(f_uni);
            std::vector<Polynomial<int64_t>> irr;
            for (size_t i = 0; i < factored.size(); ++i) {
                if (factored[i].degree() <= 0) continue;
                int exp_i = factored.exponent(i);
                for (int e = 0; e < exp_i; ++e)
                    irr.push_back(factored[i]);
            }

            bestPoints = points;
            bestUni = f_uni;
            bestFactors = irr;
            found = true;
            return;
        }

        for (const auto& val : candidates) {
            if (found) return;
            points.push_back({auxVars[idx], val});
            tryPoints(idx + 1, points);
            points.pop_back();
        }
    };

    std::vector<std::pair<int, T>> points;
    tryPoints(0, points);

    if (!found || bestFactors.size() <= 1) {
        // 既約
        std::vector<MPoly> result;
        if (cont != T(1) && cont != T(0))
            result.push_back(MPoly(cont));
        result.push_back(prim);
        return result;
    }

    // 補助変数の最大次数
    int maxAuxDeg = 0;
    for (int v : auxVars)
        for (const auto& t : prim.terms())
            maxAuxDeg = std::max(maxAuxDeg, t.mono[v]);

    // 二変数 (補助変数 1 つ) の場合: Hensel lifting + 組み合わせ
    if (auxVars.size() == 1 && bestFactors.size() >= 2) {
        int nf = static_cast<int>(bestFactors.size());

        // 因子の部分集合を試す
        for (int mask = 1; mask < (1 << nf) - 1; ++mask) {
            Polynomial<int64_t> prod_g({1}), prod_h({1});
            for (int j = 0; j < nf; ++j) {
                if (mask & (1 << j)) prod_g = prod_g * bestFactors[j];
                else prod_h = prod_h * bestFactors[j];
            }
            // Hensel の規約: deg(g0) <= deg(h0)
            if (prod_g.degree() > prod_h.degree())
                std::swap(prod_g, prod_h);

            // evalVal != 0 の場合: f(x, y) → f(x, y + a) に変換して
            // evalVal=0 で Hensel lifting し、結果を y → y - a に戻す
            T evalVal = bestPoints[0].second;
            MPoly f_shifted = prim;
            if (evalVal != T(0)) {
                // y → y + a に変換 (substitute で y=y+a はできないので手動)
                auto y_coeffs = viewAsUnivariate(prim, auxVars[0]);
                f_shifted = MPoly();
                for (int j = 0; j < static_cast<int>(y_coeffs.size()); ++j) {
                    if (y_coeffs[j].isZero()) continue;
                    // y^j → (y + a)^j = Σ C(j,m) a^{j-m} y^m
                    for (int m = 0; m <= j; ++m) {
                        int64_t binom = 1;
                        for (int i = 0; i < j - m; ++i)
                            binom = binom * (j - i) / (i + 1);
                        T a_pow = T(1);
                        for (int i = 0; i < j - m; ++i) a_pow = a_pow * evalVal;
                        auto term = y_coeffs[j] * (static_cast<T>(binom) * a_pow);
                        // y^m の項として追加
                        for (const auto& t : term.terms()) {
                            auto ev = t.mono.exponents();
                            while (static_cast<int>(ev.size()) <= auxVars[0])
                                ev.push_back(0);
                            ev[auxVars[0]] = m;
                            f_shifted.addTerm(Monomial(std::move(ev)), t.coeff);
                        }
                    }
                }
            }

            MPoly g_out, h_out;
            bool ok = detail::bivarHenselLift<T>(
                f_shifted, g_out, h_out, prod_g, prod_h,
                mainVar, auxVars[0]);

            // 結果を y → y - a に戻す
            if (ok && evalVal != T(0)) {
                auto shiftBack = [&](MPoly& poly) {
                    auto yc = viewAsUnivariate(poly, auxVars[0]);
                    MPoly result;
                    for (int j = 0; j < static_cast<int>(yc.size()); ++j) {
                        if (yc[j].isZero()) continue;
                        // y^j → (y - a)^j
                        for (int m = 0; m <= j; ++m) {
                            int64_t binom = 1;
                            for (int i = 0; i < j - m; ++i)
                                binom = binom * (j - i) / (i + 1);
                            T neg_a_pow = T(1);
                            for (int i = 0; i < j - m; ++i)
                                neg_a_pow = neg_a_pow * (T(0) - evalVal);
                            auto term = yc[j] * (static_cast<T>(binom) * neg_a_pow);
                            for (const auto& t : term.terms()) {
                                auto ev = t.mono.exponents();
                                while (static_cast<int>(ev.size()) <= auxVars[0])
                                    ev.push_back(0);
                                ev[auxVars[0]] = m;
                                result.addTerm(Monomial(std::move(ev)), t.coeff);
                            }
                        }
                    }
                    poly = result;
                };
                shiftBack(g_out);
                shiftBack(h_out);
            }

            if (ok) {
                std::vector<MPoly> result;
                if (cont != T(1) && cont != T(0))
                    result.push_back(MPoly(cont));
                auto gf = factorMultivariate(g_out);
                auto hf = factorMultivariate(h_out);
                result.insert(result.end(), gf.begin(), gf.end());
                result.insert(result.end(), hf.begin(), hf.end());
                return result;
            }
        }
    }

    // フォールバック: 既約
    std::vector<MPoly> result;
    if (cont != T(1) && cont != T(0))
        result.push_back(MPoly(cont));
    result.push_back(prim);
    return result;
}

// ================================================================
// GF(p) 上の多変数多項式因数分解
// ================================================================
//
// 手順:
//   1. 補助変数に値を代入して 1 変数化
//   2. gfp_poly::factorize で GF(p)[x] 上の因数分解
//   3. Hensel lifting で多変数因子を復元 (全て mod p)
//
// 係数型 T = int64_t で、全演算を mod p で行う。

namespace detail {

/// GF(p) 上の 1 変数多項式 (vector<int64_t>) を MultivariatePolynomial に変換
template<typename T>
MultivariatePolynomial<T> univarToMultivar(const std::vector<int64_t>& poly, int var) {
    MultivariatePolynomial<T> result;
    for (int i = 0; i < static_cast<int>(poly.size()); ++i) {
        if (poly[i] == 0) continue;
        std::vector<int> exps(var + 1, 0);
        exps[var] = i;
        result.addTerm(Monomial(exps), static_cast<T>(poly[i]));
    }
    return result;
}

/// MultivariatePolynomial の全係数を mod p にする
template<typename T>
MultivariatePolynomial<T> modPoly(const MultivariatePolynomial<T>& f, int64_t p) {
    MultivariatePolynomial<T> result;
    for (const auto& term : f.terms()) {
        T c = term.coeff;
        int64_t cv = static_cast<int64_t>(c) % p;
        if (cv < 0) cv += p;
        if (cv != 0) result.addTerm(term.mono, static_cast<T>(cv));
    }
    return result;
}

} // namespace detail

/// GF(p) 上の多変数多項式を因数分解
/// @param f  多変数多項式 (係数は [0, p) の整数)
/// @param p  素数
/// @return   既約因子のリスト (GF(p) 上)
template<typename T>
std::vector<MultivariatePolynomial<T>>
factorMultivariateModP(const MultivariatePolynomial<T>& f, int64_t p) {
    using MPoly = MultivariatePolynomial<T>;

    if (f.isZero()) return {};

    // 変数の特定
    auto vars = f.variables();
    if (vars.empty()) return {f};

    int mainVar = vars[0];
    if (vars.size() == 1) {
        // 1 変数: gfp_poly::factorize を直接使用
        // MultivariatePolynomial → vector<int64_t> に変換
        int d = 0;
        for (const auto& t : f.terms())
            d = std::max(d, t.mono[mainVar]);
        std::vector<int64_t> coeffs(d + 1, 0);
        for (const auto& t : f.terms())
            coeffs[t.mono[mainVar]] = static_cast<int64_t>(t.coeff) % p;

        auto factors = gfp_poly::factorize(coeffs, p);
        std::vector<MPoly> result;
        for (auto& [poly, mult] : factors) {
            auto mp = detail::univarToMultivar<T>(poly, mainVar);
            for (int i = 0; i < mult; ++i) result.push_back(mp);
        }
        return result;
    }

    // 多変数: 補助変数に 0 を代入して 1 変数化し、因数分解 → Hensel lifting
    // 簡易版: 既存の factorMultivariate を mod p で実行
    // (全中間演算を mod p にする完全版は将来課題)
    MPoly fm = detail::modPoly(f, p);
    auto intFactors = factorMultivariate(fm);

    // 結果の各因子を mod p に正規化
    std::vector<MPoly> result;
    for (auto& factor : intFactors) {
        auto normed = detail::modPoly(factor, p);
        if (!normed.isZero()) result.push_back(std::move(normed));
    }
    return result;
}

// ================================================================
// Zippel 疎補間法による多変数因数分解
// ================================================================
//
// 密な Hensel lifting で中間式が膨張する問題を、
// 評価・補間ベースのアプローチで回避する。
//
// アルゴリズム (簡略版):
//   1. 素数 p を選び、補助変数にランダム値を代入して 1 変数化
//   2. GF(p)[x] で因数分解 → 因子の骨格 (次数パターン) を決定
//   3. 複数の評価点で 1 変数因数分解を繰り返す
//   4. 各因子の係数を Newton 補間 (または Lagrange 補間) で復元
//   5. 有理数復元 (CRT or rational reconstruction) で Z[x,y,...] の因子を得る
//
// 参考: Zippel "Probabilistic Algorithms for Sparse Polynomials" (1979)

namespace detail {

/// Lagrange 補間 (mod p): 点列 (x_i, y_i) から多項式を復元
inline std::vector<int64_t> lagrangeInterpolate(
    const std::vector<int64_t>& xs, const std::vector<int64_t>& ys, int64_t p)
{
    int n = static_cast<int>(xs.size());
    std::vector<int64_t> result(n, 0);

    for (int i = 0; i < n; ++i) {
        if (ys[i] == 0) continue;
        // 基底多項式 L_i(x) = Π_{j≠i} (x - x_j) / (x_i - x_j)
        std::vector<int64_t> basis = {1};
        int64_t denom = 1;
        for (int j = 0; j < n; ++j) {
            if (j == i) continue;
            // basis *= (x - x_j)
            std::vector<int64_t> newBasis(basis.size() + 1, 0);
            for (size_t k = 0; k < basis.size(); ++k) {
                newBasis[k] = gfp_poly::mod(newBasis[k] + gfp_poly::mod(-xs[j] * basis[k], p), p);
                newBasis[k + 1] = gfp_poly::mod(newBasis[k + 1] + basis[k], p);
            }
            basis = newBasis;
            denom = gfp_poly::mod(denom * gfp_poly::mod(xs[i] - xs[j], p), p);
        }
        // basis *= y_i / denom
        int64_t inv_denom = power_mod(gfp_poly::mod(denom, p), p - 2, p);
        int64_t scale = gfp_poly::mod(ys[i] * inv_denom, p);
        for (size_t k = 0; k < basis.size() && k < result.size(); ++k)
            result[k] = gfp_poly::mod(result[k] + gfp_poly::mod(scale * basis[k], p), p);
    }
    return result;
}

} // namespace detail

/// Zippel 疎補間法による因数分解 (GF(p) 上)
/// @param f       多変数多項式 (int64_t 係数)
/// @param p       素数
/// @param mainVar 主変数のインデックス
/// @return        既約因子のリスト (GF(p) 上)
template<typename T>
std::vector<MultivariatePolynomial<T>>
factorZippel(const MultivariatePolynomial<T>& f, int64_t p) {
    using MPoly = MultivariatePolynomial<T>;

    if (f.isZero()) return {};

    auto vars = f.variables();
    if (vars.empty()) return {f};

    // 1 変数の場合は直接因数分解
    if (vars.size() == 1) {
        return factorMultivariateModP(f, p);
    }

    int mainVar = vars[0];
    std::vector<int> auxVars(vars.begin() + 1, vars.end());

    // mainVar の次数
    int mainDeg = 0;
    for (const auto& t : f.terms())
        mainDeg = std::max(mainDeg, t.mono[mainVar]);

    // 評価点の数 (各補助変数の最大次数 + 余裕)
    int maxAuxDeg = 0;
    for (const auto& t : f.terms())
        for (int v : auxVars)
            maxAuxDeg = std::max(maxAuxDeg, t.mono[v]);
    int numPoints = maxAuxDeg + 5;  // 余裕を持たせる

    // Step 1: 複数の評価点で 1 変数因数分解
    std::mt19937_64 rng(12345);
    struct PointResult {
        std::vector<int64_t> evalPoint;  // 各補助変数の値
        std::vector<gfp_poly::Poly> factors;
    };
    std::vector<PointResult> pointResults;

    int targetNumFactors = -1;  // 最初の評価で因子数を決定

    for (int pt = 0; pt < numPoints && pt < 50; ++pt) {
        // ランダム評価点
        std::vector<int64_t> evalPt;
        for (size_t i = 0; i < auxVars.size(); ++i)
            evalPt.push_back(1 + static_cast<int64_t>(rng() % static_cast<uint64_t>(p - 1)));

        // f(x, a1, a2, ...) を計算 → 1 変数多項式
        gfp_poly::Poly univar(mainDeg + 1, 0);
        for (const auto& t : f.terms()) {
            int64_t coeff = gfp_poly::mod(static_cast<int64_t>(t.coeff), p);
            // 補助変数の評価: Π a_v^{e_v}
            for (size_t vi = 0; vi < auxVars.size(); ++vi) {
                int exp = t.mono[auxVars[vi]];
                if (exp > 0) coeff = gfp_poly::mod(coeff * detail::power_mod(evalPt[vi], exp, p), p);
            }
            int mainExp = t.mono[mainVar];
            univar[mainExp] = gfp_poly::mod(univar[mainExp] + coeff, p);
        }
        univar = gfp_poly::normalize(univar, p);

        // 因数分解
        auto factors = gfp_poly::factorize(univar, p);

        // 重複度 1 のもののみ (square-free 部分)
        std::vector<gfp_poly::Poly> sqfFactors;
        for (auto& [poly, mult] : factors) {
            for (int i = 0; i < mult; ++i) sqfFactors.push_back(poly);
        }

        if (targetNumFactors < 0) {
            targetNumFactors = static_cast<int>(sqfFactors.size());
        } else if (static_cast<int>(sqfFactors.size()) != targetNumFactors) {
            continue;  // 因子数が一致しない評価点はスキップ
        }

        // 次数パターンでソート (因子の対応づけ)
        std::sort(sqfFactors.begin(), sqfFactors.end(),
                  [](const auto& a, const auto& b) {
                      if (a.size() != b.size()) return a.size() < b.size();
                      return a < b;
                  });

        pointResults.push_back({evalPt, std::move(sqfFactors)});
    }

    if (pointResults.empty() || targetNumFactors <= 1) {
        // 既約 or 評価点不足
        return {detail::modPoly(f, p)};
    }

    // Step 2: 各因子の係数を補間で復元
    // 簡易版: 1 つの補助変数のみ補間 (2 変数の場合)
    if (auxVars.size() == 1) {
        std::vector<MPoly> result;

        for (int fi = 0; fi < targetNumFactors; ++fi) {
            // 各次数の主変数の係数を、補助変数の値ごとに収集
            int factorDeg = gfp_poly::deg(pointResults[0].factors[fi]);

            MPoly factor;
            for (int xd = 0; xd <= factorDeg; ++xd) {
                // (evalPt, coeff_xd) のペアを収集
                std::vector<int64_t> xs, ys;
                for (auto& pr : pointResults) {
                    if (fi >= static_cast<int>(pr.factors.size())) continue;
                    xs.push_back(pr.evalPoint[0]);
                    int64_t c = (xd < static_cast<int>(pr.factors[fi].size()))
                                ? pr.factors[fi][xd] : 0;
                    ys.push_back(gfp_poly::mod(c, p));
                }

                // Lagrange 補間 → 補助変数の多項式
                auto coeffPoly = detail::lagrangeInterpolate(xs, ys, p);

                for (int yd = 0; yd < static_cast<int>(coeffPoly.size()); ++yd) {
                    if (coeffPoly[yd] == 0) continue;
                    std::vector<int> exps(std::max(mainVar, auxVars[0]) + 1, 0);
                    exps[mainVar] = xd;
                    exps[auxVars[0]] = yd;
                    factor.addTerm(Monomial(exps), static_cast<T>(coeffPoly[yd]));
                }
            }

            if (!factor.isZero()) result.push_back(std::move(factor));
        }

        // 検証: 因子の積が f (mod p) と一致するか
        if (!result.empty()) {
            MPoly product = result[0];
            for (size_t i = 1; i < result.size(); ++i) {
                product = product * result[i];
                product = detail::modPoly(product, p);
            }
            if (detail::modPoly(product, p) == detail::modPoly(f, p)) {
                return result;
            }
        }
    }

    // フォールバック: 既存の factorMultivariateModP を使用
    return factorMultivariateModP(f, p);
}

} // namespace calx
