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

// FiniteField.hpp
// 有限体 GF(p) および拡大体 GF(p^n) のクラス
//
// 提供するクラス:
//   GFp       — 素体 GF(p) のコンテキスト (p はランタイム指定の素数)
//   GFpElement — GF(p) の元
//   GFpn      — 拡大体 GF(p^n) のコンテキスト (既約多項式で定義)
//   GFpnElement — GF(p^n) の元 (多項式表現)
//
// 用途: 符号理論、楕円曲線暗号、代数的整数論
//
// 設計:
//   - GFp/GFpn はフィールドコンテキスト (法、既約多項式を保持)
//   - Element は値 + コンテキストへのポインタ
//   - int64_t ベース (p < 2^31 程度の素数に対応)
//   - 既存の ModularInt<P> (コンパイル時法) との差別化:
//     ランタイム法指定、拡大体対応、原始根探索
//
// 参照:
//   Lidl & Niederreiter, "Finite Fields" (1997)
//   FLINT: fq, fq_nmod, fq_zech (3種の表現)

#ifndef CALX_FINITE_FIELD_HPP
#define CALX_FINITE_FIELD_HPP

#include <math/core/common.hpp>
#include <vector>
#include <cstdint>
#include <cassert>
#include <stdexcept>
#include <algorithm>
#include <iosfwd>
#include <sstream>
#include <string>
#include <cmath>
#include <random>
#include <numeric>

namespace calx {

// 前方宣言
class GFp;
class GFpElement;
class GFpn;
class GFpnElement;

// ================================================================
// ユーティリティ (内部用)
// ================================================================

namespace detail {

/// 64-bit mod (常に非負)
[[nodiscard]] inline constexpr int64_t pos_mod(int64_t x, int64_t p) {
    int64_t r = x % p;
    return r < 0 ? r + p : r;
}

/// 乗算 mod (オーバーフロー安全)
/// p < 2^31 なら a*b < 2^62 でオーバーフローしない
[[nodiscard]] inline int64_t mul_mod(int64_t a, int64_t b, int64_t m) {
    return (a % m) * (b % m) % m;
}

/// 冪剰余: base^exp mod m (繰り返し二乗法)
[[nodiscard]] inline int64_t power_mod(int64_t base, int64_t exp, int64_t m) {
    if (m == 1) return 0;
    int64_t result = 1;
    base = pos_mod(base, m);
    while (exp > 0) {
        if (exp & 1)
            result = mul_mod(result, base, m);
        base = mul_mod(base, base, m);
        exp >>= 1;
    }
    return result;
}

/// Miller-Rabin 素数判定 (簡易版)
[[nodiscard]] inline bool is_prime_small(int64_t n) {
    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;
}

} // namespace detail

// ================================================================
// GFp: 素体 GF(p)
// ================================================================

/// 素体 GF(p) のコンテキスト
class GFp {
    int64_t p_;  // 素数の法

public:
    /// コンストラクタ (p は素数でなければならない)
    explicit GFp(int64_t p) : p_(p) {
        if (p < 2 || !detail::is_prime_small(p))
            throw MathError("GFp: p must be a prime >= 2");
    }

    /// 法 (素数)
    [[nodiscard]] int64_t characteristic() const { return p_; }

    /// 体の位数 = p
    [[nodiscard]] int64_t order() const { return p_; }

    /// 元を生成
    [[nodiscard]] GFpElement element(int64_t v) const;

    /// 零元
    [[nodiscard]] GFpElement zero() const;

    /// 単位元
    [[nodiscard]] GFpElement one() const;

    /// 原始根 (最小の生成元) を探索
    [[nodiscard]] GFpElement generator() const;

    /// 全元を列挙 (p が小さい場合のみ実用的)
    [[nodiscard]] std::vector<GFpElement> elements() const;

    /// 等価比較
    [[nodiscard]] bool operator==(const GFp& rhs) const { return p_ == rhs.p_; }
    [[nodiscard]] bool operator!=(const GFp& rhs) const { return p_ != rhs.p_; }
};

// ================================================================
// GFpElement: GF(p) の元
// ================================================================

/// GF(p) の元
class GFpElement {
    int64_t value_;    // [0, p)
    const GFp* field_; // 所属する体

    friend class GFp;
    friend class GFpn;
    friend class GFpnElement;

    GFpElement(int64_t v, const GFp* f)
        : value_(detail::pos_mod(v, f->characteristic())), field_(f) {}

public:
    /// 値の取得
    [[nodiscard]] int64_t value() const { return value_; }

    /// 所属する体
    [[nodiscard]] const GFp& field() const { return *field_; }

    /// 零元か
    [[nodiscard]] bool isZero() const { return value_ == 0; }

    /// 単位元か
    [[nodiscard]] bool isOne() const { return value_ == 1; }

    // ============================================================
    // 算術演算
    // ============================================================

    [[nodiscard]] GFpElement operator+(const GFpElement& rhs) const {
        assert(field_ == rhs.field_ && "Elements must belong to the same field");
        int64_t p = field_->characteristic();
        int64_t sum = value_ + rhs.value_;
        return GFpElement(sum >= p ? sum - p : sum, field_);
    }

    [[nodiscard]] GFpElement operator-(const GFpElement& rhs) const {
        assert(field_ == rhs.field_);
        int64_t p = field_->characteristic();
        int64_t diff = value_ - rhs.value_;
        return GFpElement(diff < 0 ? diff + p : diff, field_);
    }

    [[nodiscard]] GFpElement operator*(const GFpElement& rhs) const {
        assert(field_ == rhs.field_);
        int64_t p = field_->characteristic();
        return GFpElement(detail::mul_mod(value_, rhs.value_, p), field_);
    }

    [[nodiscard]] GFpElement operator/(const GFpElement& rhs) const {
        if (rhs.isZero()) throw MathError("GFpElement: division by zero");
        return *this * rhs.inverse();
    }

    [[nodiscard]] GFpElement operator-() const {
        int64_t p = field_->characteristic();
        return GFpElement(value_ == 0 ? 0 : p - value_, field_);
    }

    // 複合代入
    GFpElement& operator+=(const GFpElement& rhs) { *this = *this + rhs; return *this; }
    GFpElement& operator-=(const GFpElement& rhs) { *this = *this - rhs; return *this; }
    GFpElement& operator*=(const GFpElement& rhs) { *this = *this * rhs; return *this; }
    GFpElement& operator/=(const GFpElement& rhs) { *this = *this / rhs; return *this; }

    // ============================================================
    // 冪乗・逆元
    // ============================================================

    /// 冪乗 (非負指数)
    [[nodiscard]] GFpElement pow(int64_t n) const {
        if (n < 0) return inverse().pow(-n);
        int64_t p = field_->characteristic();
        return GFpElement(detail::power_mod(value_, n, p), field_);
    }

    /// 乗法的逆元 (拡張ユークリッド互除法)
    [[nodiscard]] GFpElement inverse() const {
        if (isZero()) throw MathError("GFpElement: inverse of zero");
        int64_t p = field_->characteristic();
        // a^{p-2} mod p (フェルマーの小定理)
        return GFpElement(detail::power_mod(value_, p - 2, p), field_);
    }

    // ============================================================
    // 比較
    // ============================================================

    [[nodiscard]] bool operator==(const GFpElement& rhs) const {
        return value_ == rhs.value_ && field_ == rhs.field_;
    }
    [[nodiscard]] bool operator!=(const GFpElement& rhs) const {
        return !(*this == rhs);
    }

    // ============================================================
    // 出力
    // ============================================================

    friend std::ostream& operator<<(std::ostream& os, const GFpElement& e) {
        return os << e.value_;
    }
};

// GFp メンバ関数の実装
inline GFpElement GFp::element(int64_t v) const {
    return GFpElement(v, this);
}

inline GFpElement GFp::zero() const {
    return GFpElement(0, this);
}

inline GFpElement GFp::one() const {
    return GFpElement(1, this);
}

inline GFpElement GFp::generator() const {
    // 最小の原始根を探索
    // g^k ≠ 1 for k = 1, ..., p-2, g^{p-1} = 1
    // p-1 の素因数 q₁,...,qₖ に対し、g^{(p-1)/qᵢ} ≠ 1 を確認
    int64_t pm1 = p_ - 1;

    // p-1 の素因数分解
    std::vector<int64_t> factors;
    {
        int64_t n = pm1;
        for (int64_t d = 2; d * d <= n; ++d) {
            if (n % d == 0) {
                factors.push_back(d);
                while (n % d == 0) n /= d;
            }
        }
        if (n > 1) factors.push_back(n);
    }

    for (int64_t g = 2; g < p_; ++g) {
        bool is_gen = true;
        for (int64_t q : factors) {
            if (detail::power_mod(g, pm1 / q, p_) == 1) {
                is_gen = false;
                break;
            }
        }
        if (is_gen) return GFpElement(g, this);
    }

    // p = 2 の場合: 1 が唯一の非零元
    return GFpElement(1, this);
}

inline std::vector<GFpElement> GFp::elements() const {
    std::vector<GFpElement> result;
    result.reserve(static_cast<size_t>(p_));
    for (int64_t i = 0; i < p_; ++i)
        result.push_back(GFpElement(i, this));
    return result;
}

// ================================================================
// GFpn: 拡大体 GF(p^n)
// ================================================================

/// 拡大体 GF(p^n) のコンテキスト
/// 多項式表現: 元は GF(p)[x]/(f(x)) の剰余類
/// f(x) は n 次の既約多項式
class GFpn {
    GFp base_;                      // 基底素体
    int n_;                         // 拡大次数
    std::vector<int64_t> irred_;    // 既約多項式の係数 (昇べき順, 長さ n+1)

    // irred_[0] + irred_[1]*x + ... + irred_[n]*x^n
    // irred_[n] = 1 (モニック) を仮定

public:
    /// コンストラクタ
    /// @param base 基底素体 GF(p)
    /// @param irredPoly 既約多項式の係数 (昇べき順、モニック、長さ n+1)
    GFpn(const GFp& base, std::vector<int64_t> irredPoly)
        : base_(base), n_(static_cast<int>(irredPoly.size()) - 1),
          irred_(std::move(irredPoly))
    {
        if (n_ < 1)
            throw MathError("GFpn: irreducible polynomial must have degree >= 1");
        // モニック化: 主係数が 1 でなければ正規化
        int64_t p = base_.characteristic();
        int64_t lc = detail::pos_mod(irred_[n_], p);
        if (lc == 0)
            throw MathError("GFpn: leading coefficient must be nonzero");
        if (lc != 1) {
            int64_t inv_lc = detail::power_mod(lc, p - 2, p);
            for (auto& c : irred_)
                c = detail::mul_mod(c, inv_lc, p);
        }
        // 係数を正規化
        for (auto& c : irred_)
            c = detail::pos_mod(c, p);
    }

    /// GF(p^2) の既約多項式を自動選択して構築
    /// @param base 基底素体 GF(p)
    /// @param n 拡大次数
    static GFpn create(const GFp& base, int n);

    /// 基底素体
    [[nodiscard]] const GFp& baseField() const { return base_; }

    /// 標数
    [[nodiscard]] int64_t characteristic() const { return base_.characteristic(); }

    /// 拡大次数
    [[nodiscard]] int degree() const { return n_; }

    /// 体の位数 = p^n
    [[nodiscard]] int64_t order() const {
        int64_t result = 1;
        for (int i = 0; i < n_; ++i) result *= base_.characteristic();
        return result;
    }

    /// 既約多項式の係数
    [[nodiscard]] const std::vector<int64_t>& irreducible() const { return irred_; }

    /// 元を生成 (係数ベクトル、昇べき順)
    [[nodiscard]] GFpnElement element(std::vector<int64_t> coeffs) const;

    /// 零元
    [[nodiscard]] GFpnElement zero() const;

    /// 単位元
    [[nodiscard]] GFpnElement one() const;

    /// 生成元 α (= x mod f(x))
    [[nodiscard]] GFpnElement alpha() const;

    /// 等価比較
    [[nodiscard]] bool operator==(const GFpn& rhs) const {
        return base_ == rhs.base_ && irred_ == rhs.irred_;
    }
    [[nodiscard]] bool operator!=(const GFpn& rhs) const { return !(*this == rhs); }

private:
    friend class GFpnElement;

    /// 多項式を既約多項式で mod
    void reduce(std::vector<int64_t>& coeffs) const {
        int64_t p = base_.characteristic();
        // coeffs の次数が n 以上なら、既約多項式で割った余りにする
        while (static_cast<int>(coeffs.size()) > n_) {
            int64_t lc = detail::pos_mod(coeffs.back(), p);
            if (lc == 0) {
                coeffs.pop_back();
                continue;
            }
            int deg = static_cast<int>(coeffs.size()) - 1;
            // x^deg = x^{deg-n} * (-irred_[0] - ... - irred_[n-1]*x^{n-1})
            for (int i = 0; i < n_; ++i) {
                coeffs[deg - n_ + i] = detail::pos_mod(
                    coeffs[deg - n_ + i] - detail::mul_mod(lc, irred_[i], p), p);
            }
            coeffs.pop_back();
        }
        // 末尾のゼロを除去
        while (!coeffs.empty() && coeffs.back() == 0)
            coeffs.pop_back();
    }

    /// 多項式乗算 mod f(x)
    [[nodiscard]] std::vector<int64_t> mul_poly(
        const std::vector<int64_t>& a,
        const std::vector<int64_t>& b) const
    {
        if (a.empty() || b.empty()) return {};
        int64_t p = base_.characteristic();
        size_t n = a.size() + b.size() - 1;
        std::vector<int64_t> result(n, 0);
        for (size_t i = 0; i < a.size(); ++i)
            for (size_t j = 0; j < b.size(); ++j)
                result[i + j] = detail::pos_mod(
                    result[i + j] + detail::mul_mod(a[i], b[j], p), p);
        reduce(result);
        return result;
    }

    /// 多項式の拡張 GCD (ベズー係数付き)
    /// gcd = s*a + t*b (mod p の多項式環上)
    struct PolyGcdResult {
        std::vector<int64_t> gcd;
        std::vector<int64_t> s;
        std::vector<int64_t> t;
    };

    [[nodiscard]] PolyGcdResult poly_gcd_ext(
        std::vector<int64_t> a,
        std::vector<int64_t> b) const
    {
        int64_t p = base_.characteristic();

        // 多項式の正規化ヘルパー
        auto normalize = [](std::vector<int64_t>& v) {
            while (!v.empty() && v.back() == 0) v.pop_back();
        };

        // 多項式除算
        auto poly_divmod = [&](const std::vector<int64_t>& dividend,
                               const std::vector<int64_t>& divisor)
            -> std::pair<std::vector<int64_t>, std::vector<int64_t>>
        {
            if (divisor.empty())
                throw MathError("GFpn: polynomial division by zero");

            std::vector<int64_t> rem = dividend;
            int divDeg = static_cast<int>(divisor.size()) - 1;
            int remDeg = static_cast<int>(rem.size()) - 1;

            if (remDeg < divDeg) return {{}, rem};

            int64_t inv_lc = detail::power_mod(divisor.back(), p - 2, p);
            std::vector<int64_t> quot(remDeg - divDeg + 1, 0);

            for (int i = remDeg - divDeg; i >= 0; --i) {
                int64_t coeff = detail::mul_mod(
                    detail::pos_mod(rem[i + divDeg], p), inv_lc, p);
                quot[i] = coeff;
                for (int j = 0; j <= divDeg; ++j)
                    rem[i + j] = detail::pos_mod(
                        rem[i + j] - detail::mul_mod(coeff, divisor[j], p), p);
            }

            normalize(rem);
            normalize(quot);
            return {quot, rem};
        };

        // 多項式の減算
        auto poly_sub_mul = [&](const std::vector<int64_t>& a_vec,
                                const std::vector<int64_t>& q,
                                const std::vector<int64_t>& b_vec)
            -> std::vector<int64_t>
        {
            // a - q*b
            std::vector<int64_t> qb;
            if (!q.empty() && !b_vec.empty()) {
                qb.resize(q.size() + b_vec.size() - 1, 0);
                for (size_t i = 0; i < q.size(); ++i)
                    for (size_t j = 0; j < b_vec.size(); ++j)
                        qb[i + j] = detail::pos_mod(
                            qb[i + j] + detail::mul_mod(q[i], b_vec[j], p), p);
            }
            size_t n = std::max(a_vec.size(), qb.size());
            std::vector<int64_t> result(n, 0);
            for (size_t i = 0; i < a_vec.size(); ++i)
                result[i] = a_vec[i];
            for (size_t i = 0; i < qb.size(); ++i)
                result[i] = detail::pos_mod(result[i] - qb[i], p);
            normalize(result);
            return result;
        };

        normalize(a);
        normalize(b);

        // 拡張 GCD
        std::vector<int64_t> old_r = a, r = b;
        std::vector<int64_t> old_s = {1}, s = {};
        std::vector<int64_t> old_t = {}, t = {1};

        while (!r.empty()) {
            auto [q, rem] = poly_divmod(old_r, r);
            old_r = r;
            r = rem;

            auto new_s = poly_sub_mul(old_s, q, s);
            old_s = s;
            s = new_s;

            auto new_t = poly_sub_mul(old_t, q, t);
            old_t = t;
            t = new_t;
        }

        // モニック化
        if (!old_r.empty() && old_r.back() != 1) {
            int64_t inv = detail::power_mod(old_r.back(), p - 2, p);
            for (auto& c : old_r) c = detail::mul_mod(c, inv, p);
            for (auto& c : old_s) c = detail::mul_mod(c, inv, p);
            for (auto& c : old_t) c = detail::mul_mod(c, inv, p);
        }

        return {old_r, old_s, old_t};
    }
};

// ================================================================
// GFpnElement: GF(p^n) の元
// ================================================================

/// GF(p^n) の元 (多項式表現)
/// a_0 + a_1·α + a_2·α² + ... + a_{n-1}·α^{n-1}
/// ここで α は f(α) = 0 の根 (f は既約多項式)
class GFpnElement {
    std::vector<int64_t> coeffs_;  // 係数 (昇べき順, 長さ ≤ n)
    const GFpn* field_;

    friend class GFpn;

    GFpnElement(std::vector<int64_t> c, const GFpn* f) : coeffs_(std::move(c)), field_(f) {
        int64_t p = field_->characteristic();
        for (auto& v : coeffs_)
            v = detail::pos_mod(v, p);
        // 末尾のゼロを除去
        while (!coeffs_.empty() && coeffs_.back() == 0)
            coeffs_.pop_back();
    }

public:
    /// 係数ベクトル
    [[nodiscard]] const std::vector<int64_t>& coefficients() const { return coeffs_; }

    /// i 番目の係数 (範囲外は 0)
    [[nodiscard]] int64_t operator[](size_t i) const {
        return i < coeffs_.size() ? coeffs_[i] : 0;
    }

    /// 所属する体
    [[nodiscard]] const GFpn& field() const { return *field_; }

    /// 零元か
    [[nodiscard]] bool isZero() const { return coeffs_.empty(); }

    /// 単位元か
    [[nodiscard]] bool isOne() const {
        return coeffs_.size() == 1 && coeffs_[0] == 1;
    }

    // ============================================================
    // 算術演算
    // ============================================================

    [[nodiscard]] GFpnElement operator+(const GFpnElement& rhs) const {
        assert(field_ == rhs.field_);
        int64_t p = field_->characteristic();
        size_t n = std::max(coeffs_.size(), rhs.coeffs_.size());
        std::vector<int64_t> result(n, 0);
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result[i] = coeffs_[i];
        for (size_t i = 0; i < rhs.coeffs_.size(); ++i)
            result[i] = detail::pos_mod(result[i] + rhs.coeffs_[i], p);
        return GFpnElement(std::move(result), field_);
    }

    [[nodiscard]] GFpnElement operator-(const GFpnElement& rhs) const {
        assert(field_ == rhs.field_);
        int64_t p = field_->characteristic();
        size_t n = std::max(coeffs_.size(), rhs.coeffs_.size());
        std::vector<int64_t> result(n, 0);
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result[i] = coeffs_[i];
        for (size_t i = 0; i < rhs.coeffs_.size(); ++i)
            result[i] = detail::pos_mod(result[i] - rhs.coeffs_[i], p);
        return GFpnElement(std::move(result), field_);
    }

    [[nodiscard]] GFpnElement operator*(const GFpnElement& rhs) const {
        assert(field_ == rhs.field_);
        auto result = field_->mul_poly(coeffs_, rhs.coeffs_);
        return GFpnElement(std::move(result), field_);
    }

    [[nodiscard]] GFpnElement operator/(const GFpnElement& rhs) const {
        if (rhs.isZero()) throw MathError("GFpnElement: division by zero");
        return *this * rhs.inverse();
    }

    [[nodiscard]] GFpnElement operator-() const {
        int64_t p = field_->characteristic();
        std::vector<int64_t> result(coeffs_.size());
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result[i] = coeffs_[i] == 0 ? 0 : p - coeffs_[i];
        return GFpnElement(std::move(result), field_);
    }

    // 複合代入
    GFpnElement& operator+=(const GFpnElement& rhs) { *this = *this + rhs; return *this; }
    GFpnElement& operator-=(const GFpnElement& rhs) { *this = *this - rhs; return *this; }
    GFpnElement& operator*=(const GFpnElement& rhs) { *this = *this * rhs; return *this; }
    GFpnElement& operator/=(const GFpnElement& rhs) { *this = *this / rhs; return *this; }

    // ============================================================
    // 冪乗・逆元
    // ============================================================

    /// 冪乗 (非負指数)
    [[nodiscard]] GFpnElement pow(int64_t n) const {
        if (n < 0) return inverse().pow(-n);
        if (n == 0) return field_->one();

        GFpnElement result = field_->one();
        GFpnElement base = *this;
        uint64_t m = static_cast<uint64_t>(n);

        while (m > 0) {
            if (m & 1) result = result * base;
            base = base * base;
            m >>= 1;
        }
        return result;
    }

    /// 乗法的逆元
    /// GF(p)[x]/(f(x)) 上の拡張 GCD で計算
    [[nodiscard]] GFpnElement inverse() const {
        if (isZero()) throw MathError("GFpnElement: inverse of zero");

        // a(x) * s(x) + f(x) * t(x) = 1 を解く
        // → s(x) が逆元
        auto result = field_->poly_gcd_ext(coeffs_, field_->irred_);
        // gcd は 1 になるはず (f が既約なので)
        return GFpnElement(std::move(result.s), field_);
    }

    /// Frobenius 写像: x → x^p
    [[nodiscard]] GFpnElement frobenius() const {
        return pow(field_->characteristic());
    }

    /// ノルム: N(a) = a * a^p * a^{p^2} * ... * a^{p^{n-1}}
    [[nodiscard]] GFpElement norm() const {
        // N(a) = a^{(p^n - 1)/(p - 1)}
        int64_t p = field_->characteristic();
        int n = field_->degree();
        int64_t pn = field_->order();
        int64_t exp = (pn - 1) / (p - 1);
        GFpnElement result = this->pow(exp);
        // 結果は GF(p) の元 (定数多項式)
        return field_->baseField().element(result.isZero() ? 0 : result[0]);
    }

    /// トレース: Tr(a) = a + a^p + a^{p^2} + ... + a^{p^{n-1}}
    [[nodiscard]] GFpElement trace() const {
        int64_t p = field_->characteristic();
        int n = field_->degree();
        GFpnElement sum = *this;
        GFpnElement current = *this;
        for (int i = 1; i < n; ++i) {
            current = current.frobenius();
            sum = sum + current;
        }
        // 結果は GF(p) の元
        return field_->baseField().element(sum.isZero() ? 0 : sum[0]);
    }

    /// 最小多項式 (GF(p) 上の)
    [[nodiscard]] std::vector<int64_t> minimalPolynomial() const {
        // α の共役は α, α^p, α^{p^2}, ..., α^{p^{d-1}}
        // 最小多項式 = Π(x - α^{p^i})
        int64_t p = field_->characteristic();
        int n = field_->degree();

        // 共役の列挙
        std::vector<GFpnElement> conjugates;
        GFpnElement current = *this;
        for (int i = 0; i < n; ++i) {
            // 既に見た共役と一致するか確認
            bool found = false;
            for (const auto& c : conjugates) {
                if (c == current) { found = true; break; }
            }
            if (found) break;
            conjugates.push_back(current);
            current = current.frobenius();
        }

        // 最小多項式 = Π(x - conjugate_i) を GF(p^n)[x] 上で計算
        // 結果の係数は GF(p) に入る
        // ここでは GFpn 上の多項式として計算
        // poly = {1} (定数 1)
        std::vector<std::vector<int64_t>> poly = {{1}};  // 1

        for (const auto& conj : conjugates) {
            // (x - conj) を掛ける
            // 新しい多項式の次数が 1 増える
            std::vector<std::vector<int64_t>> new_poly(poly.size() + 1);
            // x * poly
            for (size_t i = 0; i < poly.size(); ++i)
                new_poly[i + 1] = poly[i];
            // - conj * poly
            for (size_t i = 0; i < poly.size(); ++i) {
                auto neg_conj_poly_i = field_->mul_poly(
                    conj.coeffs_, poly[i]);
                // 符号を反転
                for (auto& c : neg_conj_poly_i)
                    c = detail::pos_mod(-c, p);
                // 加算
                size_t len = std::max(new_poly[i].size(), neg_conj_poly_i.size());
                new_poly[i].resize(len, 0);
                for (size_t j = 0; j < neg_conj_poly_i.size(); ++j)
                    new_poly[i][j] = detail::pos_mod(
                        new_poly[i][j] + neg_conj_poly_i[j], p);
            }
            poly = new_poly;
        }

        // poly[i] は x^i の係数 (GF(p^n) の元)
        // これを GF(p) の係数に変換 (定数項のみ取り出す)
        std::vector<int64_t> result(poly.size());
        for (size_t i = 0; i < poly.size(); ++i)
            result[i] = poly[i].empty() ? 0 : poly[i][0];

        return result;
    }

    // ============================================================
    // 比較
    // ============================================================

    [[nodiscard]] bool operator==(const GFpnElement& rhs) const {
        if (field_ != rhs.field_) return false;
        size_t n = std::max(coeffs_.size(), rhs.coeffs_.size());
        for (size_t i = 0; i < n; ++i) {
            int64_t a = i < coeffs_.size() ? coeffs_[i] : 0;
            int64_t b = i < rhs.coeffs_.size() ? rhs.coeffs_[i] : 0;
            if (a != b) return false;
        }
        return true;
    }
    [[nodiscard]] bool operator!=(const GFpnElement& rhs) const {
        return !(*this == rhs);
    }

    // ============================================================
    // 出力
    // ============================================================

    /// 文字列表現 (多項式形式)
    [[nodiscard]] std::string toString(const std::string& var = "α") const {
        if (coeffs_.empty()) return "0";

        std::string s;
        bool first = true;
        for (int i = static_cast<int>(coeffs_.size()) - 1; i >= 0; --i) {
            if (coeffs_[i] == 0) continue;

            if (!first && coeffs_[i] > 0) s += "+";

            if (i == 0) {
                s += std::to_string(coeffs_[i]);
            } else if (i == 1) {
                if (coeffs_[i] == 1) s += var;
                else s += std::to_string(coeffs_[i]) + var;
            } else {
                if (coeffs_[i] == 1) s += var + "^" + std::to_string(i);
                else s += std::to_string(coeffs_[i]) + var + "^" + std::to_string(i);
            }
            first = false;
        }
        return s.empty() ? "0" : s;
    }

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

// GFpn メンバ関数の実装
inline GFpnElement GFpn::element(std::vector<int64_t> coeffs) const {
    // n 次以上の項があれば reduce
    if (static_cast<int>(coeffs.size()) > n_)
        reduce(coeffs);
    return GFpnElement(std::move(coeffs), this);
}

inline GFpnElement GFpn::zero() const {
    return GFpnElement({}, this);
}

inline GFpnElement GFpn::one() const {
    return GFpnElement({1}, this);
}

inline GFpnElement GFpn::alpha() const {
    // α = x mod f(x)
    if (n_ == 1)
        return GFpnElement({}, this);  // GF(p) 自身の場合
    return GFpnElement({0, 1}, this);
}

inline GFpn GFpn::create(const GFp& base, int n) {
    // n 次の既約多項式を探索
    if (n < 1)
        throw MathError("GFpn::create: degree must be >= 1");

    int64_t p = base.characteristic();

    if (n == 1) {
        // GF(p) 自身: f(x) = x
        return GFpn(base, {0, 1});
    }

    // 小さい体の場合、全探索で既約多項式を見つける
    // f(x) = x^n + a_{n-1}*x^{n-1} + ... + a_0
    // 既約性テスト: GF(p)[x] 上で gcd(x^{p^i} - x, f(x)) を計算

    // 係数空間を走査
    int64_t total = 1;
    for (int i = 0; i < n; ++i) total *= p;

    for (int64_t code = 0; code < total; ++code) {
        std::vector<int64_t> poly(n + 1);
        int64_t tmp = code;
        for (int i = 0; i < n; ++i) {
            poly[i] = tmp % p;
            tmp /= p;
        }
        poly[n] = 1; // モニック

        // 既約性テスト
        // f(x) が既約 ⟺ すべての 1 ≤ k ≤ n/2 に対して gcd(x^{p^k} - x, f(x)) = 1
        // かつ x^{p^n} ≡ x (mod f(x))

        GFpn candidate(base, poly);
        bool irreducible = true;

        // x^{p^k} mod f(x) を計算して gcd を確認
        // x = {0, 1}
        std::vector<int64_t> x_pk = {0, 1};  // x

        for (int k = 1; k <= n / 2; ++k) {
            // x^{p^k} = (x^{p^{k-1}})^p mod f(x)
            // 冪乗を繰り返し二乗法で計算
            std::vector<int64_t> base_poly = x_pk;
            std::vector<int64_t> result_poly = {1};  // 1
            int64_t exp = p;
            while (exp > 0) {
                if (exp & 1)
                    result_poly = candidate.mul_poly(result_poly, base_poly);
                base_poly = candidate.mul_poly(base_poly, base_poly);
                exp >>= 1;
            }
            x_pk = result_poly;

            // gcd(x^{p^k} - x, f(x))
            std::vector<int64_t> diff = x_pk;
            // x^{p^k} - x
            if (diff.size() < 2) diff.resize(2, 0);
            diff[1] = detail::pos_mod(diff[1] - 1, p);
            // 末尾ゼロ除去
            while (!diff.empty() && diff.back() == 0) diff.pop_back();

            auto gcd_result = candidate.poly_gcd_ext(diff, poly);
            // gcd の次数が 0 (定数) でなければ可約
            if (gcd_result.gcd.size() > 1) {
                irreducible = false;
                break;
            }
        }

        if (irreducible) return candidate;
    }

    throw MathError("GFpn::create: failed to find irreducible polynomial");
}

// ================================================================
// GF(p)[x] 上の多項式因数分解
// ================================================================
//
// 3段階アルゴリズム:
//   1. SFF (Square-Free Factorization) — 重複因子の分離
//   2. DDF (Distinct-Degree Factorization) — 次数ごとに分離
//   3. EDF (Equal-Degree Factorization) — Cantor-Zassenhaus
//
// 多項式は vector<int64_t> で表現 (coeffs[i] = x^i の係数, mod p)

namespace gfp_poly {

using Poly = std::vector<int64_t>;

// --- 基本演算 (mod p) ---

inline int64_t mod(int64_t a, int64_t p) {
    a %= p; return (a < 0) ? a + p : a;
}

inline Poly normalize(Poly f, int64_t p) {
    for (auto& c : f) c = mod(c, p);
    while (f.size() > 1 && f.back() == 0) f.pop_back();
    return f;
}

inline int deg(const Poly& f) { return static_cast<int>(f.size()) - 1; }

inline Poly polyAdd(const Poly& a, const Poly& b, int64_t p) {
    Poly r(std::max(a.size(), b.size()), 0);
    for (size_t i = 0; i < a.size(); ++i) r[i] = mod(r[i] + a[i], p);
    for (size_t i = 0; i < b.size(); ++i) r[i] = mod(r[i] + b[i], p);
    return normalize(r, p);
}

inline Poly polySub(const Poly& a, const Poly& b, int64_t p) {
    Poly r(std::max(a.size(), b.size()), 0);
    for (size_t i = 0; i < a.size(); ++i) r[i] = mod(r[i] + a[i], p);
    for (size_t i = 0; i < b.size(); ++i) r[i] = mod(r[i] - b[i], p);
    return normalize(r, p);
}

inline Poly polyMul(const Poly& a, const Poly& b, int64_t p) {
    if (a.empty() || b.empty()) return {0};
    Poly r(a.size() + b.size() - 1, 0);
    for (size_t i = 0; i < a.size(); ++i)
        for (size_t j = 0; j < b.size(); ++j)
            r[i + j] = mod(r[i + j] + detail::mul_mod(a[i], b[j], p), p);
    return normalize(r, p);
}

inline Poly polyScale(const Poly& a, int64_t s, int64_t p) {
    Poly r(a.size());
    for (size_t i = 0; i < a.size(); ++i) r[i] = mod(detail::mul_mod(a[i], s, p), p);
    return normalize(r, p);
}

// 除算: a = q*b + r, return {q, r}
inline std::pair<Poly, Poly> polyDivMod(const Poly& a, const Poly& b, int64_t p) {
    if (deg(b) < 0 || (b.size() == 1 && b[0] == 0))
        throw MathError("gfp_poly::polyDivMod: division by zero");
    if (deg(a) < deg(b)) return {{0}, a};

    Poly r = a;
    int64_t inv_lc = detail::power_mod(b.back(), p - 2, p);
    int dB = deg(b);
    Poly q(r.size() - b.size() + 1, 0);

    for (int i = deg(r); i >= dB; --i) {
        if (r[i] == 0) continue;
        int64_t coeff = mod(detail::mul_mod(r[i], inv_lc, p), p);
        q[i - dB] = coeff;
        for (int j = 0; j <= dB; ++j)
            r[i - dB + j] = mod(r[i - dB + j] - detail::mul_mod(coeff, b[j], p), p);
    }
    return {normalize(q, p), normalize(r, p)};
}

inline Poly polyMod(const Poly& a, const Poly& b, int64_t p) {
    return polyDivMod(a, b, p).second;
}

inline Poly polyGcd(Poly a, Poly b, int64_t p) {
    while (deg(b) >= 0 && !(b.size() == 1 && b[0] == 0)) {
        auto r = polyMod(a, b, p);
        a = b; b = r;
    }
    if (!a.empty() && a.back() != 0) {
        int64_t inv = detail::power_mod(a.back(), p - 2, p);
        a = polyScale(a, inv, p);  // モニック化
    }
    return a;
}

// x^n mod f(x) (mod p) — 反復二乗法
inline Poly polyPowMod(Poly base, int64_t exp, const Poly& modulus, int64_t p) {
    Poly result = {1};
    base = polyMod(base, modulus, p);
    while (exp > 0) {
        if (exp & 1) result = polyMod(polyMul(result, base, p), modulus, p);
        base = polyMod(polyMul(base, base, p), modulus, p);
        exp >>= 1;
    }
    return result;
}

// 形式微分 f'(x) mod p
inline Poly polyDeriv(const Poly& f, int64_t p) {
    if (f.size() <= 1) return {0};
    Poly d(f.size() - 1);
    for (size_t i = 1; i < f.size(); ++i)
        d[i - 1] = mod(detail::mul_mod(f[i], static_cast<int64_t>(i), p), p);
    return normalize(d, p);
}

// --- SFF: Square-Free Factorization ---
// f = Π f_i^i where each f_i is square-free
inline std::vector<std::pair<Poly, int>> squareFreeFactorization(Poly f, int64_t p) {
    std::vector<std::pair<Poly, int>> factors;
    f = normalize(f, p);
    if (deg(f) <= 0) return factors;

    // モニック化
    int64_t lc = f.back();
    if (lc != 1) {
        int64_t inv = detail::power_mod(lc, p - 2, p);
        f = polyScale(f, inv, p);
    }

    Poly fp = polyDeriv(f, p);
    if (deg(fp) < 0 || (fp.size() == 1 && fp[0] == 0)) {
        // f'=0: f(x) = g(x^p) — p 乗根を取る
        Poly g(deg(f) / static_cast<int>(p) + 1, 0);
        for (int i = 0; i <= deg(f); i += static_cast<int>(p))
            g[i / static_cast<int>(p)] = f[i];
        auto sub = squareFreeFactorization(g, p);
        for (auto& [poly, mult] : sub)
            factors.push_back({poly, mult * static_cast<int>(p)});
        return factors;
    }

    Poly c = polyGcd(f, fp, p);
    Poly w = polyDivMod(f, c, p).first;

    int i = 1;
    while (deg(w) > 0) {
        Poly y = polyGcd(w, c, p);
        Poly z = polyDivMod(w, y, p).first;
        if (deg(z) > 0) factors.push_back({z, i});
        w = y;
        c = polyDivMod(c, y, p).first;
        ++i;
    }
    if (deg(c) > 0) {
        // c(x) = g(x^p)
        Poly g(deg(c) / static_cast<int>(p) + 1, 0);
        for (int j = 0; j <= deg(c); j += static_cast<int>(p))
            g[j / static_cast<int>(p)] = c[j];
        auto sub = squareFreeFactorization(g, p);
        for (auto& [poly, mult] : sub)
            factors.push_back({poly, mult * static_cast<int>(p)});
    }
    return factors;
}

// --- DDF: Distinct-Degree Factorization ---
// 入力: square-free なモニック多項式 f
// 出力: (g_d, d) のリスト — g_d は次数 d の既約因子の積
inline std::vector<std::pair<Poly, int>> distinctDegreeFactorization(Poly f, int64_t p) {
    std::vector<std::pair<Poly, int>> factors;
    Poly x = {0, 1};  // x
    Poly h = x;        // h = x^{p^i} mod f

    int d = 1;
    Poly fStar = f;

    while (2 * d <= deg(fStar)) {
        // h = h^p mod fStar
        h = polyPowMod(h, p, fStar, p);
        // g = gcd(h - x, fStar)
        Poly hx = polySub(h, x, p);
        Poly g = polyGcd(hx, fStar, p);

        if (deg(g) > 0) {
            factors.push_back({g, d});
            fStar = polyDivMod(fStar, g, p).first;
            h = polyMod(h, fStar, p);
        }
        ++d;
    }
    if (deg(fStar) > 0) {
        factors.push_back({fStar, deg(fStar)});
    }
    return factors;
}

// --- EDF: Equal-Degree Factorization (Cantor-Zassenhaus) ---
// 入力: f は次数 d の既約因子のみからなる square-free 多項式
// 出力: 既約因子のリスト
inline std::vector<Poly> equalDegreeFactorization(const Poly& f, int d, int64_t p) {
    int n = deg(f);
    if (n == d) return {f};
    if (n == 0) return {};

    // ランダム多項式を選んで分離を試みる
    std::mt19937_64 rng(42 + n * 1000 + d);
    int maxTries = 100;

    for (int trial = 0; trial < maxTries; ++trial) {
        // ランダム多項式 r ∈ GF(p)[x], deg < n
        Poly r(n);
        for (int i = 0; i < n; ++i) r[i] = static_cast<int64_t>(rng() % static_cast<uint64_t>(p));
        r = normalize(r, p);
        if (deg(r) <= 0) continue;

        Poly g;
        if (p == 2) {
            // GF(2): g = r + r^2 + r^{2^2} + ... + r^{2^{d-1}} mod f
            Poly trace = r;
            Poly ri = r;
            for (int i = 1; i < d; ++i) {
                ri = polyPowMod(ri, 2, f, p);
                trace = polyAdd(trace, ri, p);
            }
            g = polyGcd(trace, f, p);
        } else {
            // GF(p), p 奇数: g = gcd(r^{(p^d - 1)/2} - 1, f)
            int64_t exp = 1;
            for (int i = 0; i < d; ++i) exp *= p;
            exp = (exp - 1) / 2;
            Poly rp = polyPowMod(r, exp, f, p);
            Poly rp1 = polySub(rp, {1}, p);
            g = polyGcd(rp1, f, p);
        }

        if (deg(g) > 0 && deg(g) < n) {
            // 分離成功 — 再帰
            auto left = equalDegreeFactorization(g, d, p);
            auto right = equalDegreeFactorization(polyDivMod(f, g, p).first, d, p);
            left.insert(left.end(), right.begin(), right.end());
            return left;
        }
    }

    // フォールバック: f は既約
    return {f};
}

// --- 統合: GF(p)[x] 上の完全因数分解 ---
/// @param f 多項式 (coeffs[i] = x^i の係数)
/// @param p 素数
/// @return  (既約因子, 重複度) のリスト
inline std::vector<std::pair<Poly, int>> factorize(const Poly& f, int64_t p) {
    std::vector<std::pair<Poly, int>> result;
    Poly fn = normalize(Poly(f), p);
    if (deg(fn) <= 0) return result;

    // モニック化
    int64_t lc = fn.back();
    if (lc != 1 && lc != 0) {
        int64_t inv = detail::power_mod(lc, p - 2, p);
        fn = polyScale(fn, inv, p);
    }

    // SFF
    auto sqfFactors = squareFreeFactorization(fn, p);

    for (auto& [sqf, mult] : sqfFactors) {
        // DDF
        auto ddfFactors = distinctDegreeFactorization(sqf, p);

        for (auto& [ddf, d] : ddfFactors) {
            // EDF
            auto irred = equalDegreeFactorization(ddf, d, p);
            for (auto& factor : irred) {
                result.push_back({factor, mult});
            }
        }
    }
    return result;
}

} // namespace gfp_poly

} // namespace calx

#endif // CALX_FINITE_FIELD_HPP
