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

// NumberField.hpp
//
// 数体 Q(α) の演算
//
// 既約多項式 f(x) の根 α を添加した代数体 Q(α) = Q[x]/(f(x)) の元を、
// α の多項式 (次数 < deg(f)) として表現し、四則演算を行う。
//
// 主な API:
//   NumberField(minpoly)                — 数体の定義
//   NumberFieldElement(nf, coeffs)       — 元の構築
//   operator+, -, *, /                  — 四則演算
//   norm(x), trace(x)                   — ノルム、トレース
//   minimalPolynomial(x)                — 元の最小多項式
//   NumberField::generator()            — α (原始元)
//   NumberField::one(), zero()          — 単位元、零元

#ifndef CALX_CORE_NUMBER_FIELD_HPP
#define CALX_CORE_NUMBER_FIELD_HPP

#include <vector>
#include <cstdint>
#include <stdexcept>
#include <cmath>
#include <algorithm>
#include <numeric>
#include <complex>

namespace calx {

class NumberFieldElement;

/**
 * @brief 数体 Q(α) = Q[x]/(f(x))
 *
 * f(x) は既約多項式 (有理係数、モニック)
 */
class NumberField {
public:
    /// f(x) の係数 (coeffs[i] = x^i の係数) で数体を定義
    explicit NumberField(const std::vector<double>& minpoly)
        : minpoly_(minpoly) {
        if (minpoly.size() < 2)
            throw std::invalid_argument("NumberField: polynomial degree must be >= 1");
        // モニック化
        double lc = minpoly_.back();
        if (std::abs(lc) < 1e-15)
            throw std::invalid_argument("NumberField: leading coefficient is zero");
        for (auto& c : minpoly_) c /= lc;
        deg_ = static_cast<int>(minpoly_.size()) - 1;
    }

    int degree() const { return deg_; }
    const std::vector<double>& polynomial() const { return minpoly_; }

    /// 零元
    NumberFieldElement zero() const;
    /// 単位元
    NumberFieldElement one() const;
    /// 原始元 α
    NumberFieldElement generator() const;

private:
    std::vector<double> minpoly_;  // 最小多項式
    int deg_;
};

/**
 * @brief Q(α) の元
 *
 * a_0 + a_1·α + a_2·α² + ... + a_{n-1}·α^{n-1}
 * として表現 (n = deg(f))
 */
class NumberFieldElement {
public:
    NumberFieldElement(const NumberField* nf, const std::vector<double>& coeffs)
        : nf_(nf), coeffs_(coeffs) {
        coeffs_.resize(static_cast<std::size_t>(nf->degree()), 0.0);
    }

    const NumberField* field() const { return nf_; }
    const std::vector<double>& coefficients() const { return coeffs_; }
    int degree() const { return nf_->degree(); }

    bool isZero() const {
        return std::all_of(coeffs_.begin(), coeffs_.end(),
                          [](double c) { return std::abs(c) < 1e-15; });
    }

    /// double への変換 (次数 1 の場合のみ正確)
    double toDouble() const {
        // α の近似値を使って評価
        // α は minpoly の根
        double alpha = approxRoot();
        double result = 0;
        double alphaPow = 1.0;
        for (auto c : coeffs_) {
            result += c * alphaPow;
            alphaPow *= alpha;
        }
        return result;
    }

    // 四則演算
    friend NumberFieldElement operator+(const NumberFieldElement& a,
                                        const NumberFieldElement& b) {
        checkSameField(a, b);
        std::vector<double> result(static_cast<std::size_t>(a.degree()));
        for (int i = 0; i < a.degree(); ++i)
            result[static_cast<std::size_t>(i)] =
                a.coeffs_[static_cast<std::size_t>(i)] +
                b.coeffs_[static_cast<std::size_t>(i)];
        return NumberFieldElement(a.nf_, result);
    }

    friend NumberFieldElement operator-(const NumberFieldElement& a,
                                        const NumberFieldElement& b) {
        checkSameField(a, b);
        std::vector<double> result(static_cast<std::size_t>(a.degree()));
        for (int i = 0; i < a.degree(); ++i)
            result[static_cast<std::size_t>(i)] =
                a.coeffs_[static_cast<std::size_t>(i)] -
                b.coeffs_[static_cast<std::size_t>(i)];
        return NumberFieldElement(a.nf_, result);
    }

    NumberFieldElement operator-() const {
        std::vector<double> result(coeffs_.size());
        for (std::size_t i = 0; i < coeffs_.size(); ++i)
            result[i] = -coeffs_[i];
        return NumberFieldElement(nf_, result);
    }

    friend NumberFieldElement operator*(const NumberFieldElement& a,
                                        const NumberFieldElement& b) {
        checkSameField(a, b);
        int n = a.degree();
        // 多項式乗算後に f(x) で mod を取る
        std::vector<double> product(static_cast<std::size_t>(2 * n - 1), 0.0);
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
                product[static_cast<std::size_t>(i + j)] +=
                    a.coeffs_[static_cast<std::size_t>(i)] *
                    b.coeffs_[static_cast<std::size_t>(j)];

        // f(x) で mod を取る (上位から)
        return NumberFieldElement(a.nf_, reduce(product, a.nf_->polynomial()));
    }

    /// スカラー倍
    friend NumberFieldElement operator*(double c, const NumberFieldElement& a) {
        std::vector<double> result(a.coeffs_.size());
        for (std::size_t i = 0; i < a.coeffs_.size(); ++i)
            result[i] = c * a.coeffs_[i];
        return NumberFieldElement(a.nf_, result);
    }

    friend NumberFieldElement operator*(const NumberFieldElement& a, double c) {
        return c * a;
    }

    /// 逆元 (拡張ユークリッドの互除法)
    NumberFieldElement inverse() const {
        if (isZero())
            throw std::domain_error("NumberFieldElement: inverse of zero");

        // 拡張 GCD: gcd(coeffs_, minpoly) = 1 → coeffs_ * s + minpoly * t = 1
        // s が逆元
        auto [g, s, t] = polyExtGCD(coeffs_, nf_->polynomial());

        // g は定数 (GCD = 1 のスカラー倍)
        double scale = 1.0 / g[0];
        std::vector<double> result(s.size());
        for (std::size_t i = 0; i < s.size(); ++i)
            result[i] = s[i] * scale;

        result.resize(static_cast<std::size_t>(degree()), 0.0);
        return NumberFieldElement(nf_, result);
    }

    friend NumberFieldElement operator/(const NumberFieldElement& a,
                                        const NumberFieldElement& b) {
        return a * b.inverse();
    }

    friend bool operator==(const NumberFieldElement& a,
                           const NumberFieldElement& b) {
        checkSameField(a, b);
        for (int i = 0; i < a.degree(); ++i) {
            if (std::abs(a.coeffs_[static_cast<std::size_t>(i)] -
                        b.coeffs_[static_cast<std::size_t>(i)]) > 1e-10)
                return false;
        }
        return true;
    }

    friend bool operator!=(const NumberFieldElement& a,
                           const NumberFieldElement& b) {
        return !(a == b);
    }

private:
    const NumberField* nf_;
    std::vector<double> coeffs_;  // α^i の係数

    /// minpoly の実根の近似値
    double approxRoot() const {
        auto& p = nf_->polynomial();
        // Newton 法
        double x = 1.0;
        for (int iter = 0; iter < 200; ++iter) {
            double fx = 0, fpx = 0;
            double xpow = 1.0;
            for (std::size_t i = 0; i < p.size(); ++i) {
                fx += p[i] * xpow;
                if (i > 0) fpx += p[i] * static_cast<double>(i) * xpow / x;
                xpow *= x;
            }
            if (std::abs(fpx) < 1e-15) { x += 0.1; continue; }
            double dx = fx / fpx;
            x -= dx;
            if (std::abs(dx) < 1e-14) break;
        }
        return x;
    }

    static void checkSameField(const NumberFieldElement& a,
                                const NumberFieldElement& b) {
        if (a.nf_ != b.nf_)
            throw std::invalid_argument("NumberFieldElement: different fields");
    }

    /// 多項式 a を f で割った余り (次数 < deg(f))
    static std::vector<double> reduce(std::vector<double> a,
                                       const std::vector<double>& f) {
        int degF = static_cast<int>(f.size()) - 1;
        while (static_cast<int>(a.size()) > degF) {
            double lc = a.back();
            int shift = static_cast<int>(a.size()) - static_cast<int>(f.size());
            for (int i = 0; i < static_cast<int>(f.size()); ++i) {
                a[static_cast<std::size_t>(shift + i)] -= lc * f[static_cast<std::size_t>(i)];
            }
            // 上位の (ほぼゼロの) 項を除去
            while (!a.empty() && std::abs(a.back()) < 1e-14) a.pop_back();
        }
        a.resize(static_cast<std::size_t>(degF), 0.0);
        return a;
    }

    /// 多項式の拡張 GCD
    struct PolyGCDResult {
        std::vector<double> g, s, t;
    };

    static PolyGCDResult polyExtGCD(const std::vector<double>& a,
                                     const std::vector<double>& b) {
        std::vector<double> old_r = a, r = b;
        std::vector<double> old_s = {1.0}, s = {};
        std::vector<double> old_t = {}, t = {1.0};

        while (!r.empty() && !(r.size() == 1 && std::abs(r[0]) < 1e-15)) {
            auto [q, rem] = polyDivMod(old_r, r);
            old_r = r; r = rem;

            auto new_s = polySub2(old_s, polyMul2(q, s));
            old_s = s; s = new_s;

            auto new_t = polySub2(old_t, polyMul2(q, t));
            old_t = t; t = new_t;
        }

        return {old_r, old_s, old_t};
    }

    static std::pair<std::vector<double>, std::vector<double>>
    polyDivMod(std::vector<double> a, const std::vector<double>& b) {
        if (b.empty() || (b.size() == 1 && std::abs(b[0]) < 1e-15))
            throw std::invalid_argument("polyDivMod: division by zero");

        if (a.size() < b.size()) return {{}, a};

        int n = static_cast<int>(a.size()) - static_cast<int>(b.size());
        std::vector<double> q(static_cast<std::size_t>(n + 1), 0.0);

        for (int i = n; i >= 0; --i) {
            q[static_cast<std::size_t>(i)] = a[static_cast<std::size_t>(i + static_cast<int>(b.size()) - 1)] / b.back();
            for (int j = 0; j < static_cast<int>(b.size()); ++j) {
                a[static_cast<std::size_t>(i + j)] -= q[static_cast<std::size_t>(i)] * b[static_cast<std::size_t>(j)];
            }
        }

        // 余りをトリム
        a.resize(b.size() - 1);
        while (!a.empty() && std::abs(a.back()) < 1e-14) a.pop_back();

        return {q, a};
    }

    static std::vector<double> polyMul2(const std::vector<double>& a,
                                         const std::vector<double>& b) {
        if (a.empty() || b.empty()) return {};
        std::size_t n = a.size() + b.size() - 1;
        std::vector<double> result(n, 0.0);
        for (std::size_t i = 0; i < a.size(); ++i)
            for (std::size_t j = 0; j < b.size(); ++j)
                result[i + j] += a[i] * b[j];
        return result;
    }

    static std::vector<double> polySub2(const std::vector<double>& a,
                                         const std::vector<double>& b) {
        std::size_t n = std::max(a.size(), b.size());
        std::vector<double> result(n, 0.0);
        for (std::size_t i = 0; i < a.size(); ++i) result[i] += a[i];
        for (std::size_t i = 0; i < b.size(); ++i) result[i] -= b[i];
        while (!result.empty() && std::abs(result.back()) < 1e-14) result.pop_back();
        return result;
    }
};

// NumberField のメンバ関数実装
inline NumberFieldElement NumberField::zero() const {
    return NumberFieldElement(this, std::vector<double>(static_cast<std::size_t>(deg_), 0.0));
}

inline NumberFieldElement NumberField::one() const {
    std::vector<double> c(static_cast<std::size_t>(deg_), 0.0);
    c[0] = 1.0;
    return NumberFieldElement(this, c);
}

inline NumberFieldElement NumberField::generator() const {
    std::vector<double> c(static_cast<std::size_t>(deg_), 0.0);
    if (deg_ >= 2) c[1] = 1.0;
    else c[0] = 1.0;
    return NumberFieldElement(this, c);
}

// =====================================================================
// ノルムとトレース
// =====================================================================

/**
 * @brief 数体の元のトレース Tr(x)
 *
 * x の最小多項式 (または特性多項式) の根の和。
 * Q(α) の元 x = Σ a_i α^i のトレースは、
 * x の Q 上の乗法による n×n 行列のトレース。
 */
[[nodiscard]]
inline double trace(const NumberFieldElement& x) {
    int n = x.degree();
    auto& f = x.field()->polynomial();
    auto& c = x.coefficients();

    // x の乗法行列 M_x を構築 (α^i → x·α^i mod f)
    // M_x[j][i] = (x · α^i の α^j 係数)
    // トレースは M_x の対角和

    // 各 α^i の表現を計算
    // α^0, α^1, ..., α^{n-1} は基底
    // α^n = -f[0] - f[1]α - ... - f[n-1]α^{n-1}

    // α^k (k >= n) の基底表現を保持
    std::vector<std::vector<double>> alphaPowers(static_cast<std::size_t>(2 * n),
        std::vector<double>(static_cast<std::size_t>(n), 0.0));
    for (int i = 0; i < n; ++i)
        alphaPowers[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)] = 1.0;

    // α^n 以降を計算
    for (int k = n; k < 2 * n; ++k) {
        // α^k = α · α^{k-1}
        auto& prev = alphaPowers[static_cast<std::size_t>(k - 1)];
        auto& cur = alphaPowers[static_cast<std::size_t>(k)];
        // α · (Σ a_i α^i) = Σ a_i α^{i+1}
        // α^{n-1} の係数が非ゼロなら α^n の置換が必要
        for (int i = 0; i < n - 1; ++i) {
            cur[static_cast<std::size_t>(i + 1)] += prev[static_cast<std::size_t>(i)];
        }
        // α^{n-1} · α = α^n = -Σ f[i] α^i
        double topCoeff = prev[static_cast<std::size_t>(n - 1)];
        for (int i = 0; i < n; ++i) {
            cur[static_cast<std::size_t>(i)] -= topCoeff * f[static_cast<std::size_t>(i)];
        }
    }

    // x · α^i = (Σ c[j] α^j) · α^i = Σ c[j] α^{i+j}
    // トレース = Σ_i (x·α^i の α^i 係数)
    double tr = 0;
    for (int i = 0; i < n; ++i) {
        // x · α^i の α^i 係数
        double diag = 0;
        for (int j = 0; j < n; ++j) {
            int power = i + j;
            diag += c[static_cast<std::size_t>(j)] *
                    alphaPowers[static_cast<std::size_t>(power)][static_cast<std::size_t>(i)];
        }
        tr += diag;
    }
    return tr;
}

/**
 * @brief 数体の元のノルム N(x)
 *
 * x の乗法行列の行列式。
 */
[[nodiscard]]
inline double norm(const NumberFieldElement& x) {
    int n = x.degree();
    auto& f = x.field()->polynomial();
    auto& c = x.coefficients();

    // α^k の基底表現
    std::vector<std::vector<double>> alphaPowers(static_cast<std::size_t>(2 * n),
        std::vector<double>(static_cast<std::size_t>(n), 0.0));
    for (int i = 0; i < n; ++i)
        alphaPowers[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)] = 1.0;
    for (int k = n; k < 2 * n; ++k) {
        auto& prev = alphaPowers[static_cast<std::size_t>(k - 1)];
        auto& cur = alphaPowers[static_cast<std::size_t>(k)];
        for (int i = 0; i < n - 1; ++i)
            cur[static_cast<std::size_t>(i + 1)] += prev[static_cast<std::size_t>(i)];
        double topCoeff = prev[static_cast<std::size_t>(n - 1)];
        for (int i = 0; i < n; ++i)
            cur[static_cast<std::size_t>(i)] -= topCoeff * f[static_cast<std::size_t>(i)];
    }

    // 乗法行列 M[row][col] = (x · α^col の α^row 係数)
    std::vector<std::vector<double>> M(static_cast<std::size_t>(n),
        std::vector<double>(static_cast<std::size_t>(n), 0.0));
    for (int col = 0; col < n; ++col) {
        for (int j = 0; j < n; ++j) {
            int power = col + j;
            for (int row = 0; row < n; ++row) {
                M[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)] +=
                    c[static_cast<std::size_t>(j)] *
                    alphaPowers[static_cast<std::size_t>(power)][static_cast<std::size_t>(row)];
            }
        }
    }

    // LU 分解で行列式
    double det = 1.0;
    for (int i = 0; i < n; ++i) {
        int maxRow = i;
        double maxVal = std::abs(M[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)]);
        for (int k = i + 1; k < n; ++k) {
            if (std::abs(M[static_cast<std::size_t>(k)][static_cast<std::size_t>(i)]) > maxVal) {
                maxVal = std::abs(M[static_cast<std::size_t>(k)][static_cast<std::size_t>(i)]);
                maxRow = k;
            }
        }
        if (maxVal < 1e-15) return 0.0;
        if (maxRow != i) {
            std::swap(M[static_cast<std::size_t>(i)], M[static_cast<std::size_t>(maxRow)]);
            det = -det;
        }
        det *= M[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)];
        for (int k = i + 1; k < n; ++k) {
            double factor = M[static_cast<std::size_t>(k)][static_cast<std::size_t>(i)] /
                           M[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)];
            for (int j = i; j < n; ++j)
                M[static_cast<std::size_t>(k)][static_cast<std::size_t>(j)] -=
                    factor * M[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
        }
    }
    return det;
}

} // namespace calx

#endif // CALX_CORE_NUMBER_FIELD_HPP
