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

// MultivariatePolynomial.hpp
// 多変数多項式テンプレートクラス (疎表現)
//
// 内部表現: (単項式, 係数) ペアのソート済みベクトル
// 単項式: 指数ベクトル (例: x^2*y^3 → {2, 3})
// 設計:
//   - T = double, float, int, Int, Float, Rational 等で使用可能
//   - 疎表現 (非零項のみ保持)
//   - 3種の単項式順序: Lex, GrLex, GRevLex
//   - 加減乗、擬除算、評価、微分
//   - content, primitivePart (整数係数)

#pragma once

#include <vector>
#include <map>
#include <numeric>
#include <iosfwd>
#include <sstream>
#include <string>
#include <concepts>
#include <algorithm>
#include <initializer_list>
#include <cassert>
#include <functional>
#include <utility>
#include <cstddef>

// Polynomial.hpp の PolynomialCoeff concept を使用
#include "math/core/Polynomial.hpp"

namespace calx {

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

enum class MonomialOrder {
    Lex,      // 辞書式 (x > y > z)
    GrLex,    // 次数付き辞書式
    GRevLex   // 次数付き逆辞書式 (グレブナー基底のデフォルト)
};

// ================================================================
// Monomial (指数ベクトル)
// ================================================================

class Monomial {
private:
    std::vector<int> exponents_;  // exponents_[i] = 変数 i の指数

    // 末尾のゼロを除去
    void normalize() {
        while (!exponents_.empty() && exponents_.back() == 0)
            exponents_.pop_back();
    }

public:
    Monomial() = default;

    explicit Monomial(std::vector<int> exp) : exponents_(std::move(exp)) {
        normalize();
    }

    Monomial(std::initializer_list<int> il) : exponents_(il) {
        normalize();
    }

    // 単一変数の単項式: 変数 var の指数 exp
    static Monomial singleVar(int var, int exp = 1) {
        if (exp == 0) return Monomial();
        std::vector<int> e(var + 1, 0);
        e[var] = exp;
        return Monomial(std::move(e));
    }

    // アクセサ
    [[nodiscard]] int operator[](size_t i) const {
        return (i < exponents_.size()) ? exponents_[i] : 0;
    }

    [[nodiscard]] size_t numVariables() const { return exponents_.size(); }

    [[nodiscard]] const std::vector<int>& exponents() const { return exponents_; }

    // 総次数
    [[nodiscard]] int totalDegree() const {
        int sum = 0;
        for (int e : exponents_) sum += e;
        return sum;
    }

    // 定数単項式か (全指数 0)
    [[nodiscard]] bool isConstant() const { return exponents_.empty(); }

    // 単項式の積
    [[nodiscard]] Monomial operator*(const Monomial& rhs) const {
        size_t n = std::max(exponents_.size(), rhs.exponents_.size());
        std::vector<int> result(n, 0);
        for (size_t i = 0; i < exponents_.size(); ++i) result[i] += exponents_[i];
        for (size_t i = 0; i < rhs.exponents_.size(); ++i) result[i] += rhs.exponents_[i];
        return Monomial(std::move(result));
    }

    // 単項式の商 (指数の差、負になる場合は割り切れない)
    // 割り切れない場合は空の optional を返すべきだが、ここでは assert で保証
    [[nodiscard]] Monomial operator/(const Monomial& rhs) const {
        size_t n = std::max(exponents_.size(), rhs.exponents_.size());
        std::vector<int> result(n, 0);
        for (size_t i = 0; i < n; ++i) {
            result[i] = (*this)[static_cast<int>(i)] - rhs[i];
            assert(result[i] >= 0 && "Monomial not divisible");
        }
        return Monomial(std::move(result));
    }

    // rhs がこの単項式を割り切るか
    [[nodiscard]] bool isDivisibleBy(const Monomial& rhs) const {
        size_t n = std::max(exponents_.size(), rhs.exponents_.size());
        for (size_t i = 0; i < n; ++i) {
            if ((*this)[static_cast<int>(i)] < rhs[i]) return false;
        }
        return true;
    }

    // 比較
    [[nodiscard]] bool operator==(const Monomial& rhs) const {
        size_t n = std::max(exponents_.size(), rhs.exponents_.size());
        for (size_t i = 0; i < n; ++i) {
            if ((*this)[static_cast<int>(i)] != rhs[i]) return false;
        }
        return true;
    }

    [[nodiscard]] bool operator!=(const Monomial& rhs) const {
        return !(*this == rhs);
    }

    // 辞書式比較 (Lex order)
    [[nodiscard]] static int compareLex(const Monomial& a, const Monomial& b) {
        size_t n = std::max(a.exponents_.size(), b.exponents_.size());
        for (size_t i = 0; i < n; ++i) {
            int ai = a[i], bi = b[i];
            if (ai != bi) return (ai > bi) ? 1 : -1;
        }
        return 0;
    }

    // 次数付き辞書式比較 (GrLex order)
    [[nodiscard]] static int compareGrLex(const Monomial& a, const Monomial& b) {
        int da = a.totalDegree(), db = b.totalDegree();
        if (da != db) return (da > db) ? 1 : -1;
        return compareLex(a, b);
    }

    // 次数付き逆辞書式比較 (GRevLex order)
    [[nodiscard]] static int compareGRevLex(const Monomial& a, const Monomial& b) {
        int da = a.totalDegree(), db = b.totalDegree();
        if (da != db) return (da > db) ? 1 : -1;
        // 逆辞書式: 最後の変数から比較し、小さい方が大きい
        size_t n = std::max(a.exponents_.size(), b.exponents_.size());
        for (size_t i = n; i > 0; --i) {
            int ai = a[i - 1], bi = b[i - 1];
            if (ai != bi) return (ai < bi) ? 1 : -1;  // 逆順
        }
        return 0;
    }

    // 指定順序での比較
    [[nodiscard]] static int compare(const Monomial& a, const Monomial& 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;
    }

    // 最小公倍単項式 (LCM)
    [[nodiscard]] static Monomial lcm(const Monomial& a, const Monomial& b) {
        size_t n = std::max(a.exponents_.size(), b.exponents_.size());
        std::vector<int> result(n, 0);
        for (size_t i = 0; i < n; ++i)
            result[i] = std::max(a[i], b[i]);
        return Monomial(std::move(result));
    }

    // 最大公約単項式 (GCD)
    [[nodiscard]] static Monomial gcd(const Monomial& a, const Monomial& b) {
        size_t n = std::max(a.exponents_.size(), b.exponents_.size());
        std::vector<int> result(n, 0);
        for (size_t i = 0; i < n; ++i)
            result[i] = std::min(a[i], b[i]);
        return Monomial(std::move(result));
    }

    // 文字列表現
    [[nodiscard]] std::string toString(const std::vector<std::string>& varNames = {}) const {
        if (isConstant()) return "1";
        std::string s;
        for (size_t i = 0; i < exponents_.size(); ++i) {
            if (exponents_[i] == 0) continue;
            if (!s.empty()) s += "*";
            std::string var = (i < varNames.size()) ? varNames[i]
                            : std::string(1, static_cast<char>('x' + i));
            s += var;
            if (exponents_[i] > 1) s += "^" + std::to_string(exponents_[i]);
        }
        return s.empty() ? "1" : s;
    }

    friend std::ostream& operator<<(std::ostream& os, const Monomial& m) {
        return os << m.toString();
    }
};

// ================================================================
// Term<T>: 係数付き単項式
// ================================================================

template<typename T>
struct Term {
    T coeff;
    Monomial mono;

    Term() : coeff(T(0)) {}
    Term(T c, Monomial m) : coeff(std::move(c)), mono(std::move(m)) {}
};

// ================================================================
// MultivariatePolynomial<T> クラス
// ================================================================

template<PolynomialCoeff T>
class MultivariatePolynomial {
private:
    // 項のリスト (単項式順序で降順ソート済み)
    std::vector<Term<T>> terms_;
    MonomialOrder order_;

    // 同類項をまとめて正規化
    void normalize() {
        if (terms_.empty()) return;

        // 順序でソート (降順)
        sortTerms();

        // 同類項をまとめる
        std::vector<Term<T>> merged;
        merged.reserve(terms_.size());
        merged.push_back(std::move(terms_[0]));

        for (size_t i = 1; i < terms_.size(); ++i) {
            if (merged.back().mono == terms_[i].mono) {
                merged.back().coeff = merged.back().coeff + terms_[i].coeff;
            } else {
                merged.push_back(std::move(terms_[i]));
            }
        }

        // ゼロ係数を除去
        terms_.clear();
        for (auto& t : merged) {
            if (!(t.coeff == T(0)))
                terms_.push_back(std::move(t));
        }
    }

    // 降順ソート
    void sortTerms() {
        MonomialOrder ord = order_;
        std::sort(terms_.begin(), terms_.end(),
            [ord](const Term<T>& a, const Term<T>& b) {
                return Monomial::compare(a.mono, b.mono, ord) > 0;
            });
    }

public:
    // ============================================================
    // コンストラクタ
    // ============================================================

    // 零多項式
    MultivariatePolynomial(MonomialOrder order = MonomialOrder::GrLex)
        : order_(order) {}

    // 定数多項式
    MultivariatePolynomial(const T& c, MonomialOrder order = MonomialOrder::GrLex)
        : order_(order) {
        if (!(c == T(0)))
            terms_.push_back(Term<T>(c, Monomial()));
    }

    // Term リストから
    MultivariatePolynomial(std::vector<Term<T>> terms,
                           MonomialOrder order = MonomialOrder::GrLex)
        : terms_(std::move(terms)), order_(order) {
        normalize();
    }

    // initializer_list<Term<T>> から
    MultivariatePolynomial(std::initializer_list<Term<T>> il,
                           MonomialOrder order = MonomialOrder::GrLex)
        : terms_(il), order_(order) {
        normalize();
    }

    // 単一の項から (係数 + 単項式)
    MultivariatePolynomial(const T& coeff, const Monomial& mono,
                           MonomialOrder order = MonomialOrder::GrLex)
        : order_(order) {
        if (!(coeff == T(0)))
            terms_.push_back(Term<T>(coeff, mono));
    }

    // 項の追加 (既存の同一単項式があれば係数を加算)
    void addTerm(const Monomial& mono, const T& coeff) {
        if (coeff == T(0)) return;
        terms_.push_back(Term<T>(coeff, mono));
        normalize();
    }

    // ============================================================
    // ファクトリ関数
    // ============================================================

    // 変数 var を返す (例: variable(0) → x, variable(1) → y)
    static MultivariatePolynomial variable(int var,
                                           MonomialOrder order = MonomialOrder::GrLex) {
        return MultivariatePolynomial(T(1), Monomial::singleVar(var), order);
    }

    // 定数多項式
    static MultivariatePolynomial constant(const T& c,
                                           MonomialOrder order = MonomialOrder::GrLex) {
        return MultivariatePolynomial(c, order);
    }

    // ============================================================
    // アクセサ
    // ============================================================

    [[nodiscard]] bool isZero() const { return terms_.empty(); }

    [[nodiscard]] size_t numTerms() const { return terms_.size(); }

    [[nodiscard]] MonomialOrder monomialOrder() const { return order_; }

    [[nodiscard]] const std::vector<Term<T>>& terms() const { return terms_; }

    // 使用変数の最大インデックス + 1
    [[nodiscard]] size_t numVariables() const {
        size_t maxVar = 0;
        for (const auto& t : terms_)
            maxVar = std::max(maxVar, t.mono.numVariables());
        return maxVar;
    }

    // 総次数 (最大の単項式の総次数)
    [[nodiscard]] int totalDegree() const {
        if (terms_.empty()) return -1;
        int maxDeg = 0;
        for (const auto& t : terms_)
            maxDeg = std::max(maxDeg, t.mono.totalDegree());
        return maxDeg;
    }

    // 特定変数の最大次数
    [[nodiscard]] int degree(int var) const {
        int maxDeg = 0;
        for (const auto& t : terms_)
            maxDeg = std::max(maxDeg, t.mono[var]);
        return terms_.empty() ? -1 : maxDeg;
    }

    // 先頭項 (最大の単項式順序を持つ項)
    [[nodiscard]] const Term<T>& leadingTerm() const {
        assert(!terms_.empty() && "Zero polynomial has no leading term");
        return terms_.front();
    }

    // 先頭単項式
    [[nodiscard]] const Monomial& leadingMonomial() const {
        return leadingTerm().mono;
    }

    // 先頭係数
    [[nodiscard]] const T& leadingCoefficient() const {
        return leadingTerm().coeff;
    }

    // 定数項
    [[nodiscard]] T constantTerm() const {
        for (const auto& t : terms_) {
            if (t.mono.isConstant()) return t.coeff;
        }
        return T(0);
    }

    // 特定の単項式の係数を取得
    [[nodiscard]] T coefficient(const Monomial& mono) const {
        for (const auto& t : terms_) {
            if (t.mono == mono) return t.coeff;
        }
        return T(0);
    }

    // ============================================================
    // 比較演算子
    // ============================================================

    [[nodiscard]] bool operator==(const MultivariatePolynomial& rhs) const {
        if (terms_.size() != rhs.terms_.size()) return false;
        for (size_t i = 0; i < terms_.size(); ++i) {
            if (terms_[i].mono != rhs.terms_[i].mono) return false;
            if (!(terms_[i].coeff == rhs.terms_[i].coeff)) return false;
        }
        return true;
    }

    [[nodiscard]] bool operator!=(const MultivariatePolynomial& rhs) const {
        return !(*this == rhs);
    }

    // ============================================================
    // 単項演算子
    // ============================================================

    [[nodiscard]] MultivariatePolynomial operator+() const { return *this; }

    [[nodiscard]] MultivariatePolynomial operator-() const {
        MultivariatePolynomial result(order_);
        result.terms_.reserve(terms_.size());
        for (const auto& t : terms_)
            result.terms_.push_back(Term<T>(-t.coeff, t.mono));
        return result;
    }

    // ============================================================
    // 加減算
    // ============================================================

    [[nodiscard]] MultivariatePolynomial operator+(const MultivariatePolynomial& rhs) const {
        MultivariatePolynomial result(order_);
        result.terms_.reserve(terms_.size() + rhs.terms_.size());

        size_t i = 0, j = 0;
        while (i < terms_.size() && j < rhs.terms_.size()) {
            int cmp = Monomial::compare(terms_[i].mono, rhs.terms_[j].mono, order_);
            if (cmp > 0) {
                result.terms_.push_back(terms_[i++]);
            } else if (cmp < 0) {
                result.terms_.push_back(rhs.terms_[j++]);
            } else {
                T c = terms_[i].coeff + rhs.terms_[j].coeff;
                if (!(c == T(0)))
                    result.terms_.push_back(Term<T>(std::move(c), terms_[i].mono));
                ++i; ++j;
            }
        }
        while (i < terms_.size()) result.terms_.push_back(terms_[i++]);
        while (j < rhs.terms_.size()) result.terms_.push_back(rhs.terms_[j++]);

        return result;
    }

    [[nodiscard]] MultivariatePolynomial operator-(const MultivariatePolynomial& rhs) const {
        return *this + (-rhs);
    }

    // ============================================================
    // 乗算
    // ============================================================

    [[nodiscard]] MultivariatePolynomial operator*(const MultivariatePolynomial& rhs) const {
        if (terms_.empty() || rhs.terms_.empty())
            return MultivariatePolynomial(order_);

        MultivariatePolynomial result(order_);
        result.terms_.reserve(terms_.size() * rhs.terms_.size());

        for (const auto& a : terms_) {
            for (const auto& b : rhs.terms_) {
                result.terms_.push_back(Term<T>(a.coeff * b.coeff, a.mono * b.mono));
            }
        }
        result.normalize();
        return result;
    }

    // ============================================================
    // スカラー演算
    // ============================================================

    [[nodiscard]] MultivariatePolynomial operator+(const T& c) const {
        return *this + MultivariatePolynomial(c, order_);
    }

    [[nodiscard]] MultivariatePolynomial operator-(const T& c) const {
        return *this - MultivariatePolynomial(c, order_);
    }

    [[nodiscard]] MultivariatePolynomial operator*(const T& c) const {
        if (c == T(0)) return MultivariatePolynomial(order_);
        MultivariatePolynomial result(order_);
        result.terms_.reserve(terms_.size());
        for (const auto& t : terms_) {
            T newCoeff = t.coeff * c;
            if (!(newCoeff == T(0)))
                result.terms_.push_back(Term<T>(std::move(newCoeff), t.mono));
        }
        return result;
    }

    [[nodiscard]] MultivariatePolynomial operator/(const T& c) const {
        assert(!(c == T(0)) && "Division by zero");
        MultivariatePolynomial result(order_);
        result.terms_.reserve(terms_.size());
        for (const auto& t : terms_) {
            T newCoeff = t.coeff / c;
            if (!(newCoeff == T(0)))
                result.terms_.push_back(Term<T>(std::move(newCoeff), t.mono));
        }
        return result;
    }

    friend MultivariatePolynomial operator+(const T& lhs, const MultivariatePolynomial& rhs) {
        return rhs + lhs;
    }
    friend MultivariatePolynomial operator-(const T& lhs, const MultivariatePolynomial& rhs) {
        return (-rhs) + lhs;
    }
    friend MultivariatePolynomial operator*(const T& lhs, const MultivariatePolynomial& rhs) {
        return rhs * lhs;
    }

    // ============================================================
    // 複合代入演算子
    // ============================================================

    MultivariatePolynomial& operator+=(const MultivariatePolynomial& rhs) {
        *this = *this + rhs; return *this;
    }
    MultivariatePolynomial& operator-=(const MultivariatePolynomial& rhs) {
        *this = *this - rhs; return *this;
    }
    MultivariatePolynomial& operator*=(const MultivariatePolynomial& rhs) {
        *this = *this * rhs; return *this;
    }
    MultivariatePolynomial& operator+=(const T& c) { *this = *this + c; return *this; }
    MultivariatePolynomial& operator-=(const T& c) { *this = *this - c; return *this; }
    MultivariatePolynomial& operator*=(const T& c) { *this = *this * c; return *this; }
    MultivariatePolynomial& operator/=(const T& c) { *this = *this / c; return *this; }

    // 多項式同士の除算・剰余 (divmod を利用)
    [[nodiscard]] MultivariatePolynomial operator/(const MultivariatePolynomial& rhs) const {
        return divmod(rhs).quotient;
    }
    [[nodiscard]] MultivariatePolynomial operator%(const MultivariatePolynomial& rhs) const {
        return divmod(rhs).remainder;
    }
    MultivariatePolynomial& operator/=(const MultivariatePolynomial& rhs) {
        *this = *this / rhs; return *this;
    }
    MultivariatePolynomial& operator%=(const MultivariatePolynomial& rhs) {
        *this = *this % rhs; return *this;
    }

    // ============================================================
    // 冪乗
    // ============================================================

    [[nodiscard]] MultivariatePolynomial pow(int n) const {
        assert(n >= 0 && "Negative exponent not supported");
        if (n == 0) return MultivariatePolynomial(T(1), order_);
        MultivariatePolynomial result(T(1), order_);
        MultivariatePolynomial base = *this;
        unsigned int m = static_cast<unsigned int>(n);
        while (m > 0) {
            if (m & 1) result *= base;
            base *= base;
            m >>= 1;
        }
        return result;
    }

    // ============================================================
    // 擬除算 (Pseudo-division)
    // ============================================================
    //
    // 多変数多項式の除算: f を f_1, ..., f_s で割る
    // f = a_1*f_1 + ... + a_s*f_s + r (r は既約)
    //
    // 単一除数の場合:
    //   lc(g)^(deg(f)-deg(g)+1) * f = q * g + r

    struct DivResult {
        MultivariatePolynomial quotient;
        MultivariatePolynomial remainder;
    };

    // 単一除数による擬除算
    [[nodiscard]] DivResult pseudoDivmod(const MultivariatePolynomial& divisor) const {
        assert(!divisor.isZero() && "Division by zero polynomial");

        MultivariatePolynomial q(order_);
        MultivariatePolynomial r(order_);
        MultivariatePolynomial p = *this;

        const Monomial& ltDiv = divisor.leadingMonomial();
        const T& lcDiv = divisor.leadingCoefficient();

        while (!p.isZero() && p.leadingMonomial().isDivisibleBy(ltDiv)) {
            Monomial quotMono = p.leadingMonomial() / ltDiv;
            T quotCoeff = p.leadingCoefficient();

            // q = q * lcDiv + quotCoeff * quotMono
            q = q * lcDiv + MultivariatePolynomial(quotCoeff, quotMono, order_);
            // p = lcDiv * p - quotCoeff * quotMono * divisor
            p = lcDiv * p - MultivariatePolynomial(quotCoeff, quotMono, order_) * divisor;
        }

        r = std::move(p);
        return DivResult{std::move(q), std::move(r)};
    }

    // 体上の係数での通常除算 (double, Float 等)
    [[nodiscard]] DivResult divmod(const MultivariatePolynomial& divisor) const {
        assert(!divisor.isZero() && "Division by zero polynomial");

        MultivariatePolynomial q(order_);
        MultivariatePolynomial p = *this;

        const Monomial& ltDiv = divisor.leadingMonomial();
        const T& lcDiv = divisor.leadingCoefficient();

        while (!p.isZero() && p.leadingMonomial().isDivisibleBy(ltDiv)) {
            T quotCoeff = p.leadingCoefficient() / lcDiv;
            Monomial quotMono = p.leadingMonomial() / ltDiv;

            MultivariatePolynomial term(quotCoeff, quotMono, order_);
            q += term;
            p -= term * divisor;
        }

        return DivResult{std::move(q), std::move(p)};
    }

    // ============================================================
    // 評価
    // ============================================================

    // 全変数に値を代入して評価
    [[nodiscard]] T evaluate(const std::vector<T>& values) const {
        T result = T(0);
        for (const auto& t : terms_) {
            T termVal = t.coeff;
            for (size_t i = 0; i < t.mono.numVariables(); ++i) {
                int exp = t.mono[i];
                if (exp > 0 && i < values.size()) {
                    // 冪乗計算
                    T base = values[i];
                    T power = T(1);
                    for (int e = 0; e < exp; ++e)
                        power = power * base;
                    termVal = termVal * power;
                }
            }
            result = result + termVal;
        }
        return result;
    }

    // operator() による評価
    [[nodiscard]] T operator()(const std::vector<T>& values) const {
        return evaluate(values);
    }

    // ============================================================
    // 微分
    // ============================================================

    // 変数 var による偏微分
    [[nodiscard]] MultivariatePolynomial derivative(int var) const {
        MultivariatePolynomial result(order_);
        for (const auto& t : terms_) {
            int exp = t.mono[var];
            if (exp == 0) continue;

            T newCoeff = t.coeff * T(exp);
            std::vector<int> newExp = t.mono.exponents();
            // 指数ベクトルを必要に応じて拡張
            if (static_cast<size_t>(var) >= newExp.size())
                newExp.resize(var + 1, 0);
            newExp[var] = exp - 1;
            result.terms_.push_back(Term<T>(std::move(newCoeff), Monomial(std::move(newExp))));
        }
        result.normalize();
        return result;
    }

    // ============================================================
    // 変数置換
    // ============================================================

    // 変数 var に定数 val を代入
    [[nodiscard]] MultivariatePolynomial substitute(int var, const T& val) const {
        MultivariatePolynomial result(order_);
        for (const auto& t : terms_) {
            int exp = t.mono[var];
            T coeff = t.coeff;
            // val^exp を計算
            T power = T(1);
            for (int e = 0; e < exp; ++e)
                power = power * val;
            coeff = coeff * power;

            // 変数 var の指数を 0 にした単項式
            std::vector<int> newExp = t.mono.exponents();
            if (static_cast<size_t>(var) < newExp.size())
                newExp[var] = 0;
            result.terms_.push_back(Term<T>(std::move(coeff), Monomial(std::move(newExp))));
        }
        result.normalize();
        return result;
    }

    // ============================================================
    // 順序変更
    // ============================================================

    [[nodiscard]] MultivariatePolynomial withOrder(MonomialOrder newOrder) const {
        MultivariatePolynomial result(terms_, newOrder);
        return result;
    }

    // ============================================================
    // 文字列変換・出力
    // ============================================================

    [[nodiscard]] std::string toString(const std::vector<std::string>& varNames = {}) const {
        if (terms_.empty()) return "0";

        std::string s;
        bool first = true;

        for (const auto& t : terms_) {
            std::ostringstream oss;
            oss << t.coeff;
            std::string coefStr = oss.str();
            std::string monoStr = t.mono.toString(varNames);

            if (!first) {
                if (!coefStr.empty() && coefStr[0] != '-')
                    s += "+";
            }

            if (t.mono.isConstant()) {
                s += coefStr;
            } else if (coefStr == "1") {
                s += monoStr;
            } else if (coefStr == "-1") {
                s += "-" + monoStr;
            } else {
                s += coefStr + "*" + monoStr;
            }
            first = false;
        }

        return s.empty() ? "0" : s;
    }

    friend std::ostream& operator<<(std::ostream& os, const MultivariatePolynomial& p) {
        return os << p.toString();
    }
};

// ================================================================
// 自由関数
// ================================================================

// 変数 (ファクトリ)
template<PolynomialCoeff T = double>
[[nodiscard]] MultivariatePolynomial<T> mvVar(int var,
                                              MonomialOrder order = MonomialOrder::GrLex) {
    return MultivariatePolynomial<T>::variable(var, order);
}

// content (整数係数): 全係数の GCD
template<PolynomialCoeff T>
[[nodiscard]] T content(const MultivariatePolynomial<T>& p) {
    static_assert(std::is_integral_v<T>,
                  "content() requires integral coefficient type");
    if (p.isZero()) return T(0);
    T g = T(0);
    for (const auto& t : p.terms()) {
        T c = t.coeff;
        if (c < T(0)) c = -c;
        if (g == T(0)) g = c;
        else if (!(c == T(0))) g = std::gcd(g, c);
    }
    return (g == T(0)) ? T(1) : g;
}

// primitivePart (整数係数): p / content(p)
template<PolynomialCoeff T>
[[nodiscard]] MultivariatePolynomial<T> primitivePart(const MultivariatePolynomial<T>& p) {
    static_assert(std::is_integral_v<T>,
                  "primitivePart() requires integral coefficient type");
    T c = content(p);
    if (c == T(0) || c == T(1)) return p;
    return p / c;
}

// ================================================================
// CAS-2: 多変数多項式 GCD (擬剰余による簡易版)
// ================================================================

/// 多変数多項式の GCD (1 変数ずつ再帰的に計算)
/// 簡易版: 辞書順で最初の変数に対して擬除算を繰り返す
template<PolynomialCoeff T>
[[nodiscard]] MultivariatePolynomial<T> multivariateGcd(
    MultivariatePolynomial<T> a,
    MultivariatePolynomial<T> b)
{
    // 零チェック
    if (a.isZero()) return b;
    if (b.isZero()) return a;

    // 定数同士
    if (a.terms().size() == 1 && b.terms().size() == 1 &&
        a.totalDegree() == 0 && b.totalDegree() == 0) {
        T g = std::gcd(a.leadCoeff(), b.leadCoeff());
        MultivariatePolynomial<T> result;
        result.addTerm(typename MultivariatePolynomial<T>::Monomial{}, g);
        return result;
    }

    // content-and-primitive-part 方式
    // (完全な実装には Brown/EEZGCD が必要 — ここでは簡易版)
    T ca = content(a), cb = content(b);
    T g = std::gcd(ca, cb);

    if (ca != T(0) && ca != T(1)) a = a / ca;
    if (cb != T(0) && cb != T(1)) b = b / cb;

    // ユークリッド互除法 (擬剰余)
    constexpr int maxIter = 100;
    for (int iter = 0; iter < maxIter; ++iter) {
        if (b.isZero()) break;
        auto r = a % b;  // operator% が定義されている場合
        a = b;
        b = r;
    }

    // content を復元
    return a * g;
}

// ================================================================
// CAS-6: 多項式の関数分解 f(x) = g(h(x))
// ================================================================

/// 1 変数多項式の関数分解 (Kozen-Landau アルゴリズムの簡易版)
/// f = g(h(x)) を探す。見つからなければ空のペアを返す。
template<PolynomialCoeff T>
[[nodiscard]] std::pair<Polynomial<T>, Polynomial<T>>
functionalDecompose(const Polynomial<T>& f) {
    int n = f.degree();
    if (n <= 1) return {};

    // n = r * s (非自明な因数) を探す
    for (int s = 2; s <= n / 2; ++s) {
        if (n % s != 0) continue;
        int r = n / s;

        // h(x) = x^s + h_{s-1}*x^{s-1} + ... を仮定
        // f(x) = g(h(x)) から g, h の係数を逐次決定
        Polynomial<T> h({T(0), T(1)});  // h = x initially
        // h = x^s + lower terms (簡易: h = x^s のみ試行)
        std::vector<T> hCoeffs(s + 1, T(0));
        hCoeffs[s] = f[n] != T(0) ? T(1) : T(0);  // 正規化
        if (hCoeffs[s] == T(0)) continue;

        // g(y) = a_r * y^r + ... を f の最高次から決定
        // g_r = f_n / h_s^r → g_r = f_n (h_s = 1)
        // 簡易版: h = x^s, g は f の s 次おきの係数
        bool valid = true;
        std::vector<T> gCoeffs(r + 1, T(0));
        for (int i = 0; i <= r; ++i) {
            int idx = i * s;
            if (idx <= n)
                gCoeffs[i] = f[idx];
        }

        // 検証: g(h(x)) == f(x) ?
        Polynomial<T> hPoly(std::vector<T>(s + 1, T(0)));
        hPoly[s] = T(1);
        Polynomial<T> gPoly(gCoeffs);
        // g(h(x)) を計算
        Polynomial<T> composed(T(0));
        Polynomial<T> hPow(T(1));
        for (int i = 0; i <= r; ++i) {
            composed = composed + hPow * gCoeffs[i];
            if (i < r) hPow = hPow * hPoly;
        }

        if (composed.degree() == f.degree()) {
            bool match = true;
            for (int i = 0; i <= n; ++i)
                if (composed[i] != f[i]) { match = false; break; }
            if (match) return {gPoly, hPoly};
        }
    }
    return {};
}

// ================================================================
// CAS-11: 置換群の基本演算
// ================================================================

/// 置換の群演算 (合成、逆元、位数)
class PermutationGroup {
public:
    using Perm = std::vector<int>;  // σ[i] = σ(i)

    /// 合成 σ ∘ τ
    static Perm compose(const Perm& sigma, const Perm& tau) {
        int n = static_cast<int>(sigma.size());
        Perm result(n);
        for (int i = 0; i < n; ++i)
            result[i] = sigma[tau[i]];
        return result;
    }

    /// 逆置換
    static Perm inverse(const Perm& sigma) {
        int n = static_cast<int>(sigma.size());
        Perm result(n);
        for (int i = 0; i < n; ++i)
            result[sigma[i]] = i;
        return result;
    }

    /// 位数 (恒等置換に戻るまでの最小合成回数)
    static int order(const Perm& sigma) {
        int n = static_cast<int>(sigma.size());
        Perm current(n);
        std::iota(current.begin(), current.end(), 0);
        for (int k = 1; k <= 1000; ++k) {
            current = compose(sigma, current);
            bool identity = true;
            for (int i = 0; i < n; ++i)
                if (current[i] != i) { identity = false; break; }
            if (identity) return k;
        }
        return -1;
    }

    /// 巡回分解
    static std::vector<std::vector<int>> cycles(const Perm& sigma) {
        int n = static_cast<int>(sigma.size());
        std::vector<bool> visited(n, false);
        std::vector<std::vector<int>> result;
        for (int i = 0; i < n; ++i) {
            if (visited[i]) continue;
            std::vector<int> cycle;
            int j = i;
            while (!visited[j]) {
                visited[j] = true;
                cycle.push_back(j);
                j = sigma[j];
            }
            if (cycle.size() > 1) result.push_back(std::move(cycle));
        }
        return result;
    }

    /// 符号 (+1 = 偶置換, -1 = 奇置換)
    static int sign(const Perm& sigma) {
        auto cyc = cycles(sigma);
        int transpositions = 0;
        for (auto& c : cyc) transpositions += static_cast<int>(c.size()) - 1;
        return (transpositions % 2 == 0) ? 1 : -1;
    }

    /// 生成元の集合から群全体を生成
    static std::vector<Perm> generate(const std::vector<Perm>& generators) {
        int n = generators.empty() ? 0 : static_cast<int>(generators[0].size());
        std::vector<Perm> group;
        Perm identity(n);
        std::iota(identity.begin(), identity.end(), 0);
        group.push_back(identity);

        bool changed = true;
        while (changed) {
            changed = false;
            std::size_t sz = group.size();
            for (std::size_t i = 0; i < sz; ++i)
                for (const auto& g : generators) {
                    auto product = compose(group[i], g);
                    bool found = false;
                    for (const auto& e : group)
                        if (e == product) { found = true; break; }
                    if (!found) { group.push_back(product); changed = true; }
                }
        }
        return group;
    }
};

// ================================================================
// CAS-3: 多変数多項式の因数分解
// ================================================================

/// 多変数多項式を 1 変数多項式として見る (指定変数以外を係数とみなす)
/// 結果: coefficients[k] = x_var^k の係数 (他の変数の多項式)
template<PolynomialCoeff T>
std::vector<MultivariatePolynomial<T>> viewAsUnivariate(
    const MultivariatePolynomial<T>& f, int var)
{
    using MPoly = MultivariatePolynomial<T>;
    int maxDeg = 0;
    for (const auto& t : f.terms())
        maxDeg = std::max(maxDeg, t.mono[var]);

    std::vector<MPoly> coeffs(maxDeg + 1);
    for (const auto& t : f.terms()) {
        int d = t.mono[var];
        auto ev = t.mono.exponents();
        if (var < static_cast<int>(ev.size())) ev[var] = 0;
        coeffs[d].addTerm(Monomial(std::move(ev)), t.coeff);
    }
    return coeffs;
}

/// 多変数多項式の変数に値を代入
template<PolynomialCoeff T>
MultivariatePolynomial<T> substitute(
    const MultivariatePolynomial<T>& f, int var, T value)
{
    using MPoly = MultivariatePolynomial<T>;
    MPoly result;
    for (const auto& t : f.terms()) {
        int d = t.mono[var];
        auto ev = t.mono.exponents();
        if (var < static_cast<int>(ev.size())) ev[var] = 0;
        T coeff = t.coeff;
        for (int i = 0; i < d; ++i) coeff *= value;
        result.addTerm(Monomial(std::move(ev)), coeff);
    }
    return result;
}

// factorMultivariate() は MultivariatePolynomial_factorization.hpp に実装

// ================================================================
// 型特性
// ================================================================

template<typename T>
struct is_multivariate_polynomial : std::false_type {};

template<PolynomialCoeff T>
struct is_multivariate_polynomial<MultivariatePolynomial<T>> : std::true_type {};

template<typename T>
inline constexpr bool is_multivariate_polynomial_v = is_multivariate_polynomial<T>::value;

} // namespace calx
