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

// GroebnerBasis.hpp — Gröbner 基底 (Buchberger アルゴリズム)
//
// 多変数多項式イデアルの Gröbner 基底を計算する。
// 単項式順序: Lex, GrLex, GrevLex をサポート。
// ヘッダオンリー、テンプレート実装。
//
// 使用例:
//   using P = calx::MvPoly<double>;
//   auto x0 = P::var(0, 2);  // x^1, 2変数
//   auto x1 = P::var(1, 2);  // y^1, 2変数
//   auto f = x0*x0 + x1*x1 - P::constant(1.0, 2);  // x^2 + y^2 - 1
//   auto g = x0 - x1;                                 // x - y
//   auto basis = calx::groebnerBasis<double>({f, g});

#ifndef CALX_MATH_CORE_GROEBNER_BASIS_HPP
#define CALX_MATH_CORE_GROEBNER_BASIS_HPP

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <functional>
#include <numeric>
#include <utility>
#include <vector>

namespace calx {

    // ================================================================
    // 単項式順序
    // ================================================================

    /// 単項式の順序 (項順序)
    enum class MonomialOrder {
        Lex,      // 辞書式順序
        GrLex,    // 次数付き辞書式順序
        GrevLex   // 次数付き逆辞書式順序
    };

    // ================================================================
    // Monomial — 単項式の指数ベクトル
    // ================================================================

    /// 単項式 x0^a0 * x1^a1 * ... を指数ベクトルで表現
    template<typename T>
    struct Monomial {
        std::vector<int> exponents;

        Monomial() = default;

        explicit Monomial(std::vector<int> exp) : exponents(std::move(exp)) {}

        explicit Monomial(size_t numVars) : exponents(numVars, 0) {}

        /// 変数の数
        size_t numVars() const { return exponents.size(); }

        /// 全次数 (指数の和)
        int totalDegree() const {
            int d = 0;
            for (int e : exponents) d += e;
            return d;
        }

        /// ゼロ単項式 (定数項) か
        bool isConstant() const {
            for (int e : exponents)
                if (e != 0) return false;
            return true;
        }

        /// 単項式の乗算 (指数の和)
        Monomial operator*(const Monomial& other) const {
            size_t n = std::max(exponents.size(), other.exponents.size());
            Monomial result(n);
            for (size_t i = 0; i < n; ++i) {
                int a = (i < exponents.size()) ? exponents[i] : 0;
                int b = (i < other.exponents.size()) ? other.exponents[i] : 0;
                result.exponents[i] = a + b;
            }
            return result;
        }

        /// 単項式の除算 (指数の差) — 割り切れる前提
        Monomial operator/(const Monomial& other) const {
            size_t n = std::max(exponents.size(), other.exponents.size());
            Monomial result(n);
            for (size_t i = 0; i < n; ++i) {
                int a = (i < exponents.size()) ? exponents[i] : 0;
                int b = (i < other.exponents.size()) ? other.exponents[i] : 0;
                result.exponents[i] = a - b;
            }
            return result;
        }

        /// other が this を割り切るか (全ての指数で other ≤ this)
        bool isDivisibleBy(const Monomial& other) const {
            size_t n = std::max(exponents.size(), other.exponents.size());
            for (size_t i = 0; i < n; ++i) {
                int a = (i < exponents.size()) ? exponents[i] : 0;
                int b = (i < other.exponents.size()) ? other.exponents[i] : 0;
                if (a < b) return false;
            }
            return true;
        }

        /// 最小公倍単項式 (各指数の max)
        static Monomial lcm(const Monomial& a, const Monomial& b) {
            size_t n = std::max(a.exponents.size(), b.exponents.size());
            Monomial result(n);
            for (size_t i = 0; i < n; ++i) {
                int ea = (i < a.exponents.size()) ? a.exponents[i] : 0;
                int eb = (i < b.exponents.size()) ? b.exponents[i] : 0;
                result.exponents[i] = std::max(ea, eb);
            }
            return result;
        }

        bool operator==(const Monomial& other) const {
            size_t n = std::max(exponents.size(), other.exponents.size());
            for (size_t i = 0; i < n; ++i) {
                int a = (i < exponents.size()) ? exponents[i] : 0;
                int b = (i < other.exponents.size()) ? other.exponents[i] : 0;
                if (a != b) return false;
            }
            return true;
        }

        bool operator!=(const Monomial& other) const { return !(*this == other); }
    };

    // ================================================================
    // 単項式比較関数
    // ================================================================

    namespace detail_groebner {

        /// Lex 順序: 左から順に指数を比較
        template<typename T>
        int compareLex(const Monomial<T>& a, const Monomial<T>& b) {
            size_t n = std::max(a.exponents.size(), b.exponents.size());
            for (size_t i = 0; i < n; ++i) {
                int ea = (i < a.exponents.size()) ? a.exponents[i] : 0;
                int eb = (i < b.exponents.size()) ? b.exponents[i] : 0;
                if (ea > eb) return 1;
                if (ea < eb) return -1;
            }
            return 0;
        }

        /// GrLex 順序: 全次数 → Lex
        template<typename T>
        int compareGrLex(const Monomial<T>& a, const Monomial<T>& b) {
            int da = a.totalDegree(), db = b.totalDegree();
            if (da > db) return 1;
            if (da < db) return -1;
            return compareLex(a, b);
        }

        /// GrevLex 順序: 全次数 → 逆辞書式 (最後の変数から逆順に比較し、小さい方が大きい)
        template<typename T>
        int compareGrevLex(const Monomial<T>& a, const Monomial<T>& b) {
            int da = a.totalDegree(), db = b.totalDegree();
            if (da > db) return 1;
            if (da < db) return -1;
            // 最後の変数から逆順に比較 (指数が小さい方が「大きい」)
            size_t n = std::max(a.exponents.size(), b.exponents.size());
            for (size_t i = n; i-- > 0; ) {
                int ea = (i < a.exponents.size()) ? a.exponents[i] : 0;
                int eb = (i < b.exponents.size()) ? b.exponents[i] : 0;
                if (ea < eb) return 1;   // 逆
                if (ea > eb) return -1;  // 逆
            }
            return 0;
        }

        /// 順序に応じた比較
        template<typename T>
        int compareMonomials(const Monomial<T>& a, const Monomial<T>& b, MonomialOrder order) {
            switch (order) {
            case MonomialOrder::Lex:     return compareLex(a, b);
            case MonomialOrder::GrLex:   return compareGrLex(a, b);
            case MonomialOrder::GrevLex: return compareGrevLex(a, b);
            }
            return 0;
        }

        /// 係数がゼロとみなせるか (浮動小数点用)
        template<typename T>
        bool isZeroCoeff(const T& c) {
            return c == T(0);
        }

        // double/float 特殊化: 許容誤差付き
        template<>
        inline bool isZeroCoeff<double>(const double& c) {
            return std::abs(c) < 1e-12;
        }

        template<>
        inline bool isZeroCoeff<float>(const float& c) {
            return std::abs(c) < 1e-6f;
        }

    } // namespace detail_groebner

    // ================================================================
    // MvPoly — 多変数多項式
    // ================================================================

    /// 多変数多項式 (係数型 T)
    /// terms はソート済みを保証 (先頭が最大項 = 先導項)
    template<typename T>
    struct MvPoly {
        std::vector<std::pair<T, Monomial<T>>> terms;  // (係数, 単項式) のリスト

        MvPoly() = default;

        /// 単項からの構築
        MvPoly(T coeff, Monomial<T> mono) {
            if (!detail_groebner::isZeroCoeff(coeff))
                terms.push_back({ std::move(coeff), std::move(mono) });
        }

        /// ゼロ多項式か
        bool isZero() const { return terms.empty(); }

        /// 定数多項式を作成
        static MvPoly constant(T val, size_t numVars) {
            return MvPoly(std::move(val), Monomial<T>(numVars));
        }

        /// 単一変数 x_i (numVars 変数中)
        static MvPoly var(size_t varIdx, size_t numVars) {
            Monomial<T> m(numVars);
            m.exponents[varIdx] = 1;
            return MvPoly(T(1), std::move(m));
        }

        /// 項を順序に従ってソート (先頭が最大)
        void sortTerms(MonomialOrder order) {
            std::sort(terms.begin(), terms.end(),
                [order](const auto& a, const auto& b) {
                    return detail_groebner::compareMonomials(a.second, b.second, order) > 0;
                });
        }

        /// 同類項をまとめ、ゼロ係数の項を除去
        void normalize(MonomialOrder order) {
            if (terms.empty()) return;
            sortTerms(order);

            std::vector<std::pair<T, Monomial<T>>> merged;
            merged.push_back(terms[0]);

            for (size_t i = 1; i < terms.size(); ++i) {
                if (terms[i].second == merged.back().second) {
                    merged.back().first = merged.back().first + terms[i].first;
                } else {
                    merged.push_back(terms[i]);
                }
            }

            terms.clear();
            for (auto& t : merged) {
                if (!detail_groebner::isZeroCoeff(t.first))
                    terms.push_back(std::move(t));
            }
        }

        /// 先導項 (最大の単項式)
        const Monomial<T>& leadingMonomial() const {
            assert(!isZero());
            return terms[0].second;
        }

        /// 先導係数
        const T& leadingCoeff() const {
            assert(!isZero());
            return terms[0].first;
        }

        /// 加算
        MvPoly operator+(const MvPoly& other) const {
            MvPoly result;
            result.terms.reserve(terms.size() + other.terms.size());
            result.terms.insert(result.terms.end(), terms.begin(), terms.end());
            result.terms.insert(result.terms.end(), other.terms.begin(), other.terms.end());
            return result;
        }

        /// 減算
        MvPoly operator-(const MvPoly& other) const {
            MvPoly result;
            result.terms.reserve(terms.size() + other.terms.size());
            result.terms.insert(result.terms.end(), terms.begin(), terms.end());
            for (auto& t : other.terms)
                result.terms.push_back({ T(0) - t.first, t.second });
            return result;
        }

        /// 乗算
        MvPoly operator*(const MvPoly& other) const {
            MvPoly result;
            result.terms.reserve(terms.size() * other.terms.size());
            for (auto& a : terms) {
                for (auto& b : other.terms) {
                    result.terms.push_back({ a.first * b.first, a.second * b.second });
                }
            }
            return result;
        }

        /// スカラー乗算
        MvPoly operator*(const T& scalar) const {
            MvPoly result;
            result.terms.reserve(terms.size());
            for (auto& t : terms)
                result.terms.push_back({ t.first * scalar, t.second });
            return result;
        }

        /// 単項式による乗算
        MvPoly multiplyByMonomial(const T& coeff, const Monomial<T>& mono) const {
            MvPoly result;
            result.terms.reserve(terms.size());
            for (auto& t : terms)
                result.terms.push_back({ t.first * coeff, t.second * mono });
            return result;
        }
    };

    // ================================================================
    // S 多項式
    // ================================================================

    /// S 多項式 S(f, g)
    /// S(f,g) = (lcm(LM(f),LM(g)) / LT(f)) * f - (lcm(LM(f),LM(g)) / LT(g)) * g
    template<typename T>
    MvPoly<T> sPolynomial(const MvPoly<T>& f, const MvPoly<T>& g,
                          MonomialOrder order = MonomialOrder::GrevLex) {
        assert(!f.isZero() && !g.isZero());

        auto lcmMono = Monomial<T>::lcm(f.leadingMonomial(), g.leadingMonomial());
        auto quotF = lcmMono / f.leadingMonomial();
        auto quotG = lcmMono / g.leadingMonomial();

        T cf = f.leadingCoeff();
        T cg = g.leadingCoeff();

        // S(f,g) = (1/cf) * quotF * f - (1/cg) * quotG * g
        auto termF = f.multiplyByMonomial(T(1) / cf, quotF);
        auto termG = g.multiplyByMonomial(T(1) / cg, quotG);

        auto result = termF - termG;
        result.normalize(order);
        return result;
    }

    // ================================================================
    // 多項式の簡約 (除算)
    // ================================================================

    /// f を多項式集合 G で簡約 (余りを返す)
    template<typename T>
    MvPoly<T> reduce(const MvPoly<T>& f, const std::vector<MvPoly<T>>& G,
                     MonomialOrder order = MonomialOrder::GrevLex) {
        MvPoly<T> remainder;
        MvPoly<T> h = f;
        h.normalize(order);

        while (!h.isZero()) {
            bool reduced = false;
            for (auto& g : G) {
                if (g.isZero()) continue;
                if (h.leadingMonomial().isDivisibleBy(g.leadingMonomial())) {
                    auto quotMono = h.leadingMonomial() / g.leadingMonomial();
                    T quotCoeff = h.leadingCoeff() / g.leadingCoeff();

                    // h = h - quotCoeff * quotMono * g
                    auto sub = g.multiplyByMonomial(quotCoeff, quotMono);
                    h = h - sub;
                    h.normalize(order);
                    reduced = true;
                    break;
                }
            }
            if (!reduced) {
                // 先導項を余りに移す
                remainder.terms.push_back(h.terms[0]);
                h.terms.erase(h.terms.begin());
                h.normalize(order);
            }
        }

        remainder.normalize(order);
        return remainder;
    }

    // ================================================================
    // Buchberger アルゴリズム
    // ================================================================

    /// Buchberger アルゴリズムで Gröbner 基底を計算
    ///
    /// 入力: 多項式イデアルの生成元
    /// 出力: 同じイデアルを生成する Gröbner 基底
    template<typename T>
    std::vector<MvPoly<T>> groebnerBasis(
        const std::vector<MvPoly<T>>& generators,
        MonomialOrder order = MonomialOrder::GrevLex)
    {
        // 入力を正規化
        std::vector<MvPoly<T>> G;
        for (auto& gen : generators) {
            MvPoly<T> p = gen;
            p.normalize(order);
            if (!p.isZero())
                G.push_back(std::move(p));
        }
        if (G.empty()) return G;

        // 未処理のペア (i, j)
        using Pair = std::pair<size_t, size_t>;
        std::vector<Pair> pairs;
        for (size_t i = 0; i < G.size(); ++i)
            for (size_t j = i + 1; j < G.size(); ++j)
                pairs.push_back({ i, j });

        while (!pairs.empty()) {
            auto [i, j] = pairs.back();
            pairs.pop_back();

            // Buchberger の第1判定基準:
            // lcm(LM(f), LM(g)) == LM(f)*LM(g) (互いに素) なら S 多項式は 0 に簡約
            auto lcmMono = Monomial<T>::lcm(G[i].leadingMonomial(), G[j].leadingMonomial());
            auto prodMono = G[i].leadingMonomial() * G[j].leadingMonomial();
            if (lcmMono == prodMono)
                continue;

            auto sp = sPolynomial(G[i], G[j], order);
            auto rem = reduce(sp, G, order);

            if (!rem.isZero()) {
                // 先導係数で正規化 (monic にする)
                T lc = rem.leadingCoeff();
                if (lc != T(1)) {
                    for (auto& t : rem.terms)
                        t.first = t.first / lc;
                }

                size_t newIdx = G.size();
                for (size_t k = 0; k < newIdx; ++k)
                    pairs.push_back({ k, newIdx });
                G.push_back(std::move(rem));
            }
        }

        // 最小化: 他の先導項で割り切れる生成元を除去
        std::vector<bool> redundant(G.size(), false);
        for (size_t i = 0; i < G.size(); ++i) {
            if (redundant[i]) continue;
            for (size_t j = 0; j < G.size(); ++j) {
                if (i == j || redundant[j]) continue;
                if (G[i].leadingMonomial().isDivisibleBy(G[j].leadingMonomial())) {
                    redundant[i] = true;
                    break;
                }
            }
        }

        std::vector<MvPoly<T>> minimal;
        for (size_t i = 0; i < G.size(); ++i) {
            if (!redundant[i])
                minimal.push_back(std::move(G[i]));
        }

        // 被約化: 各多項式を他の全多項式で簡約
        for (size_t i = 0; i < minimal.size(); ++i) {
            std::vector<MvPoly<T>> others;
            for (size_t j = 0; j < minimal.size(); ++j)
                if (j != i) others.push_back(minimal[j]);
            minimal[i] = reduce(minimal[i], others, order);
            minimal[i].normalize(order);
            // monic に正規化
            if (!minimal[i].isZero()) {
                T lc = minimal[i].leadingCoeff();
                if (lc != T(1)) {
                    for (auto& t : minimal[i].terms)
                        t.first = t.first / lc;
                }
            }
        }

        // ゼロ多項式を除去
        std::erase_if(minimal, [](const MvPoly<T>& p) { return p.isZero(); });

        // 順序でソート (先導項の大きい順)
        std::sort(minimal.begin(), minimal.end(),
            [order](const auto& a, const auto& b) {
                return detail_groebner::compareMonomials(
                    a.leadingMonomial(), b.leadingMonomial(), order) > 0;
            });

        return minimal;
    }

    // ================================================================
    // イデアル所属判定
    // ================================================================

    /// 多項式 f がイデアル (basis が Gröbner 基底) に属するか判定
    template<typename T>
    bool isInIdeal(const MvPoly<T>& f, const std::vector<MvPoly<T>>& basis,
                   MonomialOrder order = MonomialOrder::GrevLex) {
        auto rem = reduce(f, basis, order);
        return rem.isZero();
    }

} // namespace calx

#endif // CALX_MATH_CORE_GROEBNER_BASIS_HPP
