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

// AlgebraicNumber.hpp
//
// 代数的数 (Algebraic Numbers)
//
// 最小多項式 + 根の近似値による厳密表現。
// 四則演算は結果多項式 (resultant) で実装。
//
// 主な API:
//   AlgebraicNumber::fromRational(p, q)  — 有理数
//   AlgebraicNumber::fromPoly(coeffs, approx) — 多項式の根
//   AlgebraicNumber::sqrt(n)     — √n
//   operator+, -, *, /           — 四則演算
//   minimalPolynomial()          — 最小多項式
//   toDouble()                   — 近似値
//   degree()                     — 代数次数

#ifndef CALX_CORE_ALGEBRAIC_NUMBER_HPP
#define CALX_CORE_ALGEBRAIC_NUMBER_HPP

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

namespace calx {

namespace detail_alg {

/// 多項式の評価 (Horner 法)
inline double polyEval(const std::vector<double>& coeffs, double x) {
    // coeffs[i] = x^i の係数
    double result = 0;
    for (int i = static_cast<int>(coeffs.size()) - 1; i >= 0; --i) {
        result = result * x + coeffs[static_cast<std::size_t>(i)];
    }
    return result;
}

/// 多項式の乗算
inline std::vector<double> polyMul(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;
}

/// 多項式の加算
inline std::vector<double> polyAdd(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];
    return result;
}

/// 多項式の減算
inline std::vector<double> polySub(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];
    return result;
}

/// 多項式のスカラー倍
inline std::vector<double> polyScale(const std::vector<double>& a, double c) {
    std::vector<double> result(a.size());
    for (std::size_t i = 0; i < a.size(); ++i) result[i] = a[i] * c;
    return result;
}

/// 上位のゼロ係数を除去
inline void polyTrim(std::vector<double>& p) {
    while (p.size() > 1 && std::abs(p.back()) < 1e-12) p.pop_back();
}

/// 多項式の GCD (ユークリッド)
inline std::vector<double> polyGCD(std::vector<double> a, std::vector<double> b) {
    polyTrim(a);
    polyTrim(b);
    while (b.size() > 1 || (b.size() == 1 && std::abs(b[0]) > 1e-12)) {
        // a mod b
        while (a.size() >= b.size()) {
            double ratio = a.back() / b.back();
            std::size_t shift = a.size() - b.size();
            for (std::size_t i = 0; i < b.size(); ++i) {
                a[shift + i] -= ratio * b[i];
            }
            polyTrim(a);
        }
        std::swap(a, b);
    }
    // モニック化
    if (!a.empty() && std::abs(a.back()) > 1e-15) {
        double lc = a.back();
        for (auto& c : a) c /= lc;
    }
    return a;
}

/// Resultant (Sylvester 行列の行列式)
/// res_y(f(x,y), g(x,y)) where f, g are polynomials in y with coefficients that are polynomials in x
/// ここでは簡易版: 多項式 f(t), g(t) の resultant を数値的に計算
inline double resultant(const std::vector<double>& f, const std::vector<double>& g) {
    // Sylvester 行列の行列式
    int m = static_cast<int>(f.size()) - 1;  // deg(f)
    int n = static_cast<int>(g.size()) - 1;  // deg(g)
    if (m < 0 || n < 0) return 0.0;

    int size = m + n;
    if (size == 0) return 1.0;

    // Sylvester 行列を構築
    std::vector<std::vector<double>> mat(static_cast<std::size_t>(size),
                                         std::vector<double>(static_cast<std::size_t>(size), 0.0));

    // f の行 (n 行)
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j <= m; ++j) {
            mat[static_cast<std::size_t>(i)][static_cast<std::size_t>(i + j)] =
                f[static_cast<std::size_t>(m - j)];
        }
    }
    // g の行 (m 行)
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j <= n; ++j) {
            mat[static_cast<std::size_t>(n + i)][static_cast<std::size_t>(i + j)] =
                g[static_cast<std::size_t>(n - j)];
        }
    }

    // LU 分解で行列式
    double det = 1.0;
    for (int i = 0; i < size; ++i) {
        // ピボット選択
        int maxRow = i;
        double maxVal = std::abs(mat[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)]);
        for (int k = i + 1; k < size; ++k) {
            if (std::abs(mat[static_cast<std::size_t>(k)][static_cast<std::size_t>(i)]) > maxVal) {
                maxVal = std::abs(mat[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(mat[static_cast<std::size_t>(i)], mat[static_cast<std::size_t>(maxRow)]);
            det = -det;
        }

        det *= mat[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)];

        for (int k = i + 1; k < size; ++k) {
            double factor = mat[static_cast<std::size_t>(k)][static_cast<std::size_t>(i)] /
                           mat[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)];
            for (int j = i; j < size; ++j) {
                mat[static_cast<std::size_t>(k)][static_cast<std::size_t>(j)] -=
                    factor * mat[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
            }
        }
    }
    return det;
}

/// 多項式の根を Newton 法で精密化
inline double refineRoot(const std::vector<double>& coeffs, double approx,
                          int maxIter = 100) {
    // 導関数
    std::vector<double> deriv(coeffs.size() - 1);
    for (std::size_t i = 1; i < coeffs.size(); ++i)
        deriv[i - 1] = coeffs[i] * static_cast<double>(i);

    double x = approx;
    for (int i = 0; i < maxIter; ++i) {
        double fx = polyEval(coeffs, x);
        double fpx = polyEval(deriv, x);
        if (std::abs(fpx) < 1e-15) break;
        double dx = fx / fpx;
        x -= dx;
        if (std::abs(dx) < 1e-14 * std::max(1.0, std::abs(x))) break;
    }
    return x;
}

/// 合成多項式 f(y - c) の係数を求める (y = x + c への変数変換)
inline std::vector<double> polyCompose_shift(const std::vector<double>& f, double c) {
    // f(x - c): Taylor 展開
    int n = static_cast<int>(f.size());
    std::vector<double> result(static_cast<std::size_t>(n), 0.0);

    // f(y) = f(x + c) を展開
    // 二項定理で各項を展開
    for (int k = 0; k < n; ++k) {
        // f[k] * x^k → f[k] * (y-c)^k を展開
        // ここでは直接的にやる
        std::vector<double> term = {f[static_cast<std::size_t>(k)]};
        std::vector<double> factor = {-c, 1.0};  // (x - c)
        for (int j = 0; j < k; ++j) {
            term = polyMul(term, factor);
        }
        for (std::size_t i = 0; i < term.size() && i < result.size(); ++i) {
            result[i] += term[i];
        }
    }
    return result;
}

/// 合成 Resultant: α+β が根となる多項式を求める
/// f(α)=0, g(β)=0 のとき、h(α+β)=0 となる h を求める
/// h(t) = Res_y(f(t-y), g(y))
inline std::vector<double> addPoly(const std::vector<double>& f,
                                    const std::vector<double>& g,
                                    double approxSum) {
    // f(t-y) を y の多項式として展開し、g(y) との resultant を計算
    // 簡易実装: 根を列挙して補間
    // ここでは resultant ベースの実装

    int m = static_cast<int>(f.size()) - 1;
    int n = static_cast<int>(g.size()) - 1;
    int deg = m * n;

    // 数値的に: f の根 α_i と g の根 β_j を求め、
    // Π(t - (α_i + β_j)) を計算
    // ただし根を数値的に求める

    // Durand-Kerner で全根を求める
    auto findRoots = [](const std::vector<double>& poly) -> std::vector<std::complex<double>> {
        int d = static_cast<int>(poly.size()) - 1;
        if (d <= 0) return {};

        // モニック化
        double lc = poly[static_cast<std::size_t>(d)];
        std::vector<double> monic(poly.size());
        for (std::size_t i = 0; i < poly.size(); ++i)
            monic[i] = poly[i] / lc;

        // 初期点を円上に配置
        std::vector<std::complex<double>> roots(static_cast<std::size_t>(d));
        for (int i = 0; i < d; ++i) {
            double angle = 2.0 * 3.14159265358979323846 * i / d + 0.1;
            double r = 1.0 + std::abs(monic[0]);
            roots[static_cast<std::size_t>(i)] = std::complex<double>(
                r * std::cos(angle), r * std::sin(angle));
        }

        // Durand-Kerner 反復
        for (int iter = 0; iter < 1000; ++iter) {
            double maxDelta = 0;
            for (int i = 0; i < d; ++i) {
                std::complex<double> num(monic[static_cast<std::size_t>(d)]);
                for (int k = d - 1; k >= 0; --k) {
                    num = num * roots[static_cast<std::size_t>(i)] +
                          monic[static_cast<std::size_t>(k)];
                }

                std::complex<double> den(1.0);
                for (int j = 0; j < d; ++j) {
                    if (j != i) {
                        den *= (roots[static_cast<std::size_t>(i)] -
                                roots[static_cast<std::size_t>(j)]);
                    }
                }

                auto delta = num / den;
                roots[static_cast<std::size_t>(i)] -= delta;
                maxDelta = std::max(maxDelta, std::abs(delta));
            }
            if (maxDelta < 1e-14) break;
        }
        return roots;
    };

    auto rootsF = findRoots(f);
    auto rootsG = findRoots(g);

    // Π(t - (α_i + β_j)) を計算
    std::vector<double> result = {1.0};
    for (auto& a : rootsF) {
        for (auto& b : rootsG) {
            std::complex<double> sum = a + b;
            // 実数根のペアだけ使用
            if (std::abs(sum.imag()) < 1e-8) {
                double r = sum.real();
                result = polyMul(result, {-r, 1.0});
            } else {
                // 複素共役ペアは二次式 t² - 2Re(sum)t + |sum|²
                double re = sum.real();
                double im = sum.imag();
                result = polyMul(result, {re*re + im*im, -2*re, 1.0});
            }
        }
    }

    // 係数を丸める (整数に近い場合)
    for (auto& c : result) {
        double rounded = std::round(c);
        if (std::abs(c - rounded) < 1e-6) c = rounded;
    }

    polyTrim(result);
    return result;
}

/// 乗算多項式: α*β が根となる多項式
inline std::vector<double> mulPoly(const std::vector<double>& f,
                                    const std::vector<double>& g,
                                    double approxProd) {
    auto findRoots = [](const std::vector<double>& poly) -> std::vector<std::complex<double>> {
        int d = static_cast<int>(poly.size()) - 1;
        if (d <= 0) return {};
        double lc = poly[static_cast<std::size_t>(d)];
        std::vector<double> monic(poly.size());
        for (std::size_t i = 0; i < poly.size(); ++i) monic[i] = poly[i] / lc;

        std::vector<std::complex<double>> roots(static_cast<std::size_t>(d));
        for (int i = 0; i < d; ++i) {
            double angle = 2.0 * 3.14159265358979323846 * i / d + 0.1;
            double r = 1.0 + std::abs(monic[0]);
            roots[static_cast<std::size_t>(i)] = {r * std::cos(angle), r * std::sin(angle)};
        }
        for (int iter = 0; iter < 1000; ++iter) {
            double maxDelta = 0;
            for (int i = 0; i < d; ++i) {
                std::complex<double> num(monic[static_cast<std::size_t>(d)]);
                for (int k = d - 1; k >= 0; --k)
                    num = num * roots[static_cast<std::size_t>(i)] + monic[static_cast<std::size_t>(k)];
                std::complex<double> den(1.0);
                for (int j = 0; j < d; ++j)
                    if (j != i) den *= (roots[static_cast<std::size_t>(i)] - roots[static_cast<std::size_t>(j)]);
                auto delta = num / den;
                roots[static_cast<std::size_t>(i)] -= delta;
                maxDelta = std::max(maxDelta, std::abs(delta));
            }
            if (maxDelta < 1e-14) break;
        }
        return roots;
    };

    auto rootsF = findRoots(f);
    auto rootsG = findRoots(g);

    std::vector<double> result = {1.0};
    for (auto& a : rootsF) {
        for (auto& b : rootsG) {
            auto prod = a * b;
            if (std::abs(prod.imag()) < 1e-8) {
                result = polyMul(result, {-prod.real(), 1.0});
            } else {
                double re = prod.real(), im = prod.imag();
                result = polyMul(result, {re*re + im*im, -2*re, 1.0});
            }
        }
    }

    for (auto& c : result) {
        double rounded = std::round(c);
        if (std::abs(c - rounded) < 1e-6) c = rounded;
    }
    polyTrim(result);
    return result;
}

} // namespace detail_alg

/**
 * @brief 代数的数 (最小多項式 + 近似値)
 *
 * f(α) = 0 となる多項式 f と、根を特定するための近似値 approx を保持。
 * 有理数は 1 次多項式 (x - p/q) で表現。
 */
class AlgebraicNumber {
public:
    /// ゼロ
    AlgebraicNumber() : coeffs_({0.0, 1.0}), approx_(0.0) {}

    /// 有理数 p/q
    static AlgebraicNumber fromRational(int64_t p, int64_t q = 1) {
        if (q == 0) throw std::invalid_argument("AlgebraicNumber: denominator is zero");
        double val = static_cast<double>(p) / static_cast<double>(q);
        // 最小多項式: q*x - p
        return AlgebraicNumber({static_cast<double>(-p), static_cast<double>(q)}, val);
    }

    /// 多項式の根から構築
    static AlgebraicNumber fromPoly(const std::vector<double>& coeffs, double approx) {
        if (coeffs.size() < 2)
            throw std::invalid_argument("AlgebraicNumber: polynomial degree must be >= 1");
        AlgebraicNumber result(coeffs, approx);
        result.approx_ = detail_alg::refineRoot(result.coeffs_, approx);
        return result;
    }

    /// √n
    static AlgebraicNumber sqrt(int64_t n) {
        if (n < 0) throw std::domain_error("AlgebraicNumber::sqrt: negative argument");
        if (n == 0) return fromRational(0);
        // 完全平方数チェック
        int64_t s = static_cast<int64_t>(std::sqrt(static_cast<double>(n)));
        if (s * s == n) return fromRational(s);
        // 最小多項式: x² - n
        return AlgebraicNumber({static_cast<double>(-n), 0.0, 1.0}, std::sqrt(static_cast<double>(n)));
    }

    /// 代数次数
    int degree() const { return static_cast<int>(coeffs_.size()) - 1; }

    /// 最小多項式
    const std::vector<double>& minimalPolynomial() const { return coeffs_; }

    /// 近似値
    double toDouble() const { return approx_; }

    /// 符号反転
    AlgebraicNumber operator-() const {
        // f(α)=0 → f(-β)=0 where β=-α
        // g(x) = f(-x): 奇数次の係数を反転
        std::vector<double> g(coeffs_.size());
        for (std::size_t i = 0; i < coeffs_.size(); ++i) {
            g[i] = (i % 2 == 1) ? -coeffs_[i] : coeffs_[i];
        }
        return AlgebraicNumber(g, -approx_);
    }

    /// 加算
    friend AlgebraicNumber operator+(const AlgebraicNumber& a, const AlgebraicNumber& b) {
        if (a.degree() == 1 && b.degree() == 1) {
            // 有理数同士
            return fromRational(0).fromDouble(a.approx_ + b.approx_);
        }

        double approxSum = a.approx_ + b.approx_;

        if (a.degree() == 1) {
            // 有理数 + 代数的数: α が f の根 → α+r は f(x-r) の根
            double r = a.approx_;
            auto shifted = detail_alg::polyCompose_shift(b.coeffs_, r);
            return AlgebraicNumber(shifted,
                                   detail_alg::refineRoot(shifted, approxSum));
        }
        if (b.degree() == 1) {
            double r = b.approx_;
            auto shifted = detail_alg::polyCompose_shift(a.coeffs_, r);
            return AlgebraicNumber(shifted,
                                   detail_alg::refineRoot(shifted, approxSum));
        }

        // 一般的な場合: resultant
        auto h = detail_alg::addPoly(a.coeffs_, b.coeffs_, approxSum);
        // 最小多項式を求める (h の因子で approxSum を根に持つもの)
        auto minPoly = findMinimalFactor(h, approxSum);
        return AlgebraicNumber(minPoly,
                               detail_alg::refineRoot(minPoly, approxSum));
    }

    /// 減算
    friend AlgebraicNumber operator-(const AlgebraicNumber& a, const AlgebraicNumber& b) {
        return a + (-b);
    }

    /// 乗算
    friend AlgebraicNumber operator*(const AlgebraicNumber& a, const AlgebraicNumber& b) {
        if (a.degree() == 1 && b.degree() == 1) {
            return fromRational(0).fromDouble(a.approx_ * b.approx_);
        }

        double approxProd = a.approx_ * b.approx_;

        if (a.degree() == 1) {
            // 有理数 * 代数的数: f(x/r) * r^deg が結果
            double r = a.approx_;
            if (std::abs(r) < 1e-15) return AlgebraicNumber();
            std::vector<double> scaled(b.coeffs_.size());
            double rPow = 1.0;
            for (std::size_t i = 0; i < b.coeffs_.size(); ++i) {
                scaled[i] = b.coeffs_[b.coeffs_.size() - 1 - i] * rPow;
                rPow *= r;
            }
            // 反転
            std::reverse(scaled.begin(), scaled.end());
            // 実際は b(x/r) * r^deg
            scaled.clear();
            scaled.resize(b.coeffs_.size());
            double rinv = 1.0 / r;
            double rinvPow = 1.0;
            for (std::size_t i = 0; i < b.coeffs_.size(); ++i) {
                scaled[i] = b.coeffs_[i] * rinvPow;
                rinvPow *= rinv;
            }
            return AlgebraicNumber(scaled,
                                   detail_alg::refineRoot(scaled, approxProd));
        }
        if (b.degree() == 1) {
            return b * a;
        }

        auto h = detail_alg::mulPoly(a.coeffs_, b.coeffs_, approxProd);
        auto minPoly = findMinimalFactor(h, approxProd);
        return AlgebraicNumber(minPoly,
                               detail_alg::refineRoot(minPoly, approxProd));
    }

    /// 逆元
    AlgebraicNumber inverse() const {
        if (std::abs(approx_) < 1e-15)
            throw std::domain_error("AlgebraicNumber: inverse of zero");
        // f(α)=0, g(x) = x^n * f(1/x) は 1/α の最小多項式
        std::vector<double> rev(coeffs_.rbegin(), coeffs_.rend());
        return AlgebraicNumber(rev,
                               detail_alg::refineRoot(rev, 1.0 / approx_));
    }

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

    /// 等価判定 (近似値ベース)
    friend bool operator==(const AlgebraicNumber& a, const AlgebraicNumber& b) {
        return std::abs(a.approx_ - b.approx_) < 1e-10;
    }

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

    friend bool operator<(const AlgebraicNumber& a, const AlgebraicNumber& b) {
        return a.approx_ < b.approx_ - 1e-10;
    }

private:
    std::vector<double> coeffs_;  // 最小多項式 (coeffs_[i] = x^i の係数)
    double approx_;               // 近似値

    AlgebraicNumber(const std::vector<double>& coeffs, double approx)
        : coeffs_(coeffs), approx_(approx) {
        // モニック化
        if (!coeffs_.empty() && std::abs(coeffs_.back()) > 1e-15) {
            double lc = coeffs_.back();
            for (auto& c : coeffs_) c /= lc;
        }
    }

    /// double 値から有理数として構築 (内部用)
    static AlgebraicNumber fromDouble(double v) {
        return AlgebraicNumber({-v, 1.0}, v);
    }

    /// 多項式 h の因子のうち、approx を根に持つ最小次のものを求める
    static std::vector<double> findMinimalFactor(const std::vector<double>& h, double approx) {
        // 簡易版: h 自体が既約の場合はそのまま返す
        // 次数 1 なら確認して返す
        if (h.size() <= 2) return h;

        // Newton 法で根を精密化
        double root = detail_alg::refineRoot(h, approx);

        // h を (x - root) で割って、最小多項式を探す
        // 因子分解を試みる
        std::vector<double> current = h;
        detail_alg::polyTrim(current);

        // 各次数で割れるか試す
        for (int targetDeg = 1; targetDeg < static_cast<int>(current.size()) - 1; ++targetDeg) {
            // targetDeg 次の因子を探す (root が根)
            // 実際には GCD ベースの方法で最小多項式を求める

            // (x - root) で割れるか
            if (targetDeg == 1) {
                std::vector<double> lin = {-root, 1.0};
                // h を (x - root) で割る
                int n = static_cast<int>(current.size()) - 1;
                std::vector<double> quotient(static_cast<std::size_t>(n), 0.0);
                std::vector<double> rem = current;
                for (int i = n; i >= 1; --i) {
                    quotient[static_cast<std::size_t>(i - 1)] = rem[static_cast<std::size_t>(i)];
                    rem[static_cast<std::size_t>(i - 1)] += root * rem[static_cast<std::size_t>(i)];
                    rem[static_cast<std::size_t>(i)] = 0;
                }
                if (std::abs(rem[0]) < 1e-8) {
                    return lin;
                }
            }
        }

        // 最小多項式が見つからなければ h をそのまま返す
        return current;
    }
};

/**
 * @brief 代数的数の絶対値
 */
[[nodiscard]]
inline double abs(const AlgebraicNumber& x) {
    return std::abs(x.toDouble());
}

} // namespace calx

#endif // CALX_CORE_ALGEBRAIC_NUMBER_HPP
