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

// RationalFunction.hpp
// 有理関数テンプレートクラス (多項式の商 P(x)/Q(x))
// lib++20 有理式.h からの移植 (2026-02-21)
//
// 設計:
//   - 分子・分母を Polynomial<T> で保持
//   - T = double, float, Int, Float, Rational, Complex<T> 等で使用可能
//   - GCD による自動約分
//   - Horner 法による O(n) 評価

#pragma once

#include <math/core/Polynomial.hpp>

namespace calx {

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

template<PolynomialCoeff T>
class RationalFunction {
private:
    Polynomial<T> numerator_;    // 分子
    Polynomial<T> denominator_;  // 分母

public:
    using value_type = T;

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

    // デフォルト: 0/1
    RationalFunction() : numerator_(T(0)), denominator_(T(1)) {}

    // スカラーから: c/1
    RationalFunction(const T& c) : numerator_(c), denominator_(T(1)) {}

    // 多項式から: P(x)/1
    RationalFunction(const Polynomial<T>& p)
        : numerator_(p), denominator_(T(1)) {}

    RationalFunction(Polynomial<T>&& p)
        : numerator_(std::move(p)), denominator_(T(1)) {}

    // 分子・分母を指定
    RationalFunction(const Polynomial<T>& num, const Polynomial<T>& den,
                     bool doSimplify = true)
        : numerator_(num), denominator_(den)
    {
        normalizeSign_();
        if (doSimplify) simplify();
    }

    RationalFunction(Polynomial<T>&& num, Polynomial<T>&& den,
                     bool doSimplify = true)
        : numerator_(std::move(num)), denominator_(std::move(den))
    {
        normalizeSign_();
        if (doSimplify) simplify();
    }

    // 初期化子リストから (昇べきの順 — Polynomial と同じ)
    // 例: RationalFunction<double>({1,2,3}, {0,1}) → (1+2x+3x²)/(x)
    RationalFunction(std::initializer_list<T> num,
                     std::initializer_list<T> den,
                     bool doSimplify = true)
        : numerator_(num), denominator_(den)
    {
        normalizeSign_();
        if (doSimplify) simplify();
    }

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

    [[nodiscard]] const Polynomial<T>& numerator() const { return numerator_; }
    [[nodiscard]] const Polynomial<T>& denominator() const { return denominator_; }

    // 分子の次数
    [[nodiscard]] int numeratorDegree() const { return numerator_.degree(); }

    // 分母の次数
    [[nodiscard]] int denominatorDegree() const { return denominator_.degree(); }

    // 零か (分子が零)
    [[nodiscard]] bool isZero() const { return numerator_.isZero(); }

    // ============================================================
    // 約分 (GCD で分子・分母を簡約化)
    // ============================================================

    bool simplify() {
        if (numerator_.isZero()) {
            denominator_ = Polynomial<T>(T(1));
            return true;
        }

        Polynomial<T> g = gcd(numerator_, denominator_);
        if (g.degree() >= 1 || !(g.constantTerm() == T(1))) {
            numerator_ = numerator_ / g;
            denominator_ = denominator_ / g;
            return true;
        }
        return false;
    }

    /// 近似約分 (ε-ユークリッド法による近似 GCD で分子・分母を簡約化)
    ///
    /// 浮動小数点係数の有理関数で、丸め誤差により厳密 GCD が 1 に
    /// なってしまう場合に使用する。IIR フィルタの最小実現 (minimal
    /// realization) — 極零相殺による次数低減に有用。
    ///
    /// @param tolerance 剰余の係数最大絶対値がこれ以下なら零とみなす
    /// @return 約分が行われた場合 true
    bool approximateSimplify(double tolerance = 1e-10) {
        if (numerator_.isZero()) {
            denominator_ = Polynomial<T>(T(1));
            return true;
        }

        Polynomial<T> g = approximateGCD(numerator_, denominator_, tolerance);
        if (g.degree() >= 1) {
            numerator_ = numerator_ / g;
            denominator_ = denominator_ / g;
            // 除算の丸め誤差で微小な高次項が残る場合があるので再正規化
            normalizeSign_();
            return true;
        }
        return false;
    }

    // ============================================================
    // 評価 (Horner 法)
    // ============================================================

    template<typename U = T>
    [[nodiscard]] U operator()(const U& x) const {
        return numerator_(x) / denominator_(x);
    }

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

    [[nodiscard]] bool operator==(const RationalFunction& rhs) const {
        // P1/Q1 == P2/Q2 ⟺ P1*Q2 == P2*Q1
        return numerator_ * rhs.denominator_ == rhs.numerator_ * denominator_;
    }

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

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

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

    [[nodiscard]] RationalFunction operator-() const {
        return RationalFunction(-numerator_, denominator_, false);
    }

    // ============================================================
    // RationalFunction + RationalFunction
    // P1/Q1 + P2/Q2 = (P1*Q2 + P2*Q1) / (Q1*Q2)
    // ============================================================

    [[nodiscard]] RationalFunction operator+(const RationalFunction& rhs) const {
        if (denominator_ == rhs.denominator_) {
            return RationalFunction(numerator_ + rhs.numerator_, denominator_);
        }
        return RationalFunction(
            numerator_ * rhs.denominator_ + rhs.numerator_ * denominator_,
            denominator_ * rhs.denominator_
        );
    }

    [[nodiscard]] RationalFunction operator-(const RationalFunction& rhs) const {
        if (denominator_ == rhs.denominator_) {
            return RationalFunction(numerator_ - rhs.numerator_, denominator_);
        }
        return RationalFunction(
            numerator_ * rhs.denominator_ - rhs.numerator_ * denominator_,
            denominator_ * rhs.denominator_
        );
    }

    // P1/Q1 * P2/Q2 = (P1*P2) / (Q1*Q2)
    [[nodiscard]] RationalFunction operator*(const RationalFunction& rhs) const {
        return RationalFunction(
            numerator_ * rhs.numerator_,
            denominator_ * rhs.denominator_
        );
    }

    // P1/Q1 / P2/Q2 = (P1*Q2) / (Q1*P2)
    [[nodiscard]] RationalFunction operator/(const RationalFunction& rhs) const {
        return RationalFunction(
            numerator_ * rhs.denominator_,
            denominator_ * rhs.numerator_
        );
    }

    // ============================================================
    // RationalFunction + スカラー
    // ============================================================

    [[nodiscard]] RationalFunction operator+(const T& c) const {
        return RationalFunction(numerator_ + c * denominator_, denominator_);
    }
    [[nodiscard]] RationalFunction operator-(const T& c) const {
        return RationalFunction(numerator_ - c * denominator_, denominator_);
    }
    [[nodiscard]] RationalFunction operator*(const T& c) const {
        return RationalFunction(numerator_ * c, denominator_);
    }
    [[nodiscard]] RationalFunction operator/(const T& c) const {
        return RationalFunction(numerator_, denominator_ * c);
    }

    friend RationalFunction operator+(const T& c, const RationalFunction& f) {
        return f + c;
    }
    friend RationalFunction operator-(const T& c, const RationalFunction& f) {
        return RationalFunction(c * f.denominator_ - f.numerator_, f.denominator_);
    }
    friend RationalFunction operator*(const T& c, const RationalFunction& f) {
        return f * c;
    }
    friend RationalFunction operator/(const T& c, const RationalFunction& f) {
        return RationalFunction(c * f.denominator_, f.numerator_);
    }

    // ============================================================
    // RationalFunction + Polynomial
    // ============================================================

    [[nodiscard]] RationalFunction operator+(const Polynomial<T>& p) const {
        return RationalFunction(numerator_ + p * denominator_, denominator_);
    }
    [[nodiscard]] RationalFunction operator-(const Polynomial<T>& p) const {
        return RationalFunction(numerator_ - p * denominator_, denominator_);
    }
    [[nodiscard]] RationalFunction operator*(const Polynomial<T>& p) const {
        return RationalFunction(numerator_ * p, denominator_);
    }
    [[nodiscard]] RationalFunction operator/(const Polynomial<T>& p) const {
        return RationalFunction(numerator_, denominator_ * p);
    }

    friend RationalFunction operator+(const Polynomial<T>& p, const RationalFunction& f) {
        return f + p;
    }
    friend RationalFunction operator-(const Polynomial<T>& p, const RationalFunction& f) {
        return RationalFunction(p * f.denominator_ - f.numerator_, f.denominator_);
    }
    friend RationalFunction operator*(const Polynomial<T>& p, const RationalFunction& f) {
        return f * p;
    }
    friend RationalFunction operator/(const Polynomial<T>& p, const RationalFunction& f) {
        return RationalFunction(p * f.denominator_, f.numerator_);
    }

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

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

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

    RationalFunction& operator+=(const Polynomial<T>& p) { return *this = *this + p; }
    RationalFunction& operator-=(const Polynomial<T>& p) { return *this = *this - p; }
    RationalFunction& operator*=(const Polynomial<T>& p) { return *this = *this * p; }
    RationalFunction& operator/=(const Polynomial<T>& p) { return *this = *this / p; }

    // ============================================================
    // べき乗
    // ============================================================

    [[nodiscard]] RationalFunction pow(int n) const {
        if (n == 0) return RationalFunction(T(1));
        if (n > 0) {
            return RationalFunction(numerator_.pow(n), denominator_.pow(n), false);
        }
        // n < 0: 逆数のべき乗
        return RationalFunction(denominator_.pow(-n), numerator_.pow(-n), false);
    }

    // ============================================================
    // 文字列変換・ストリーム出力
    // ============================================================

    [[nodiscard]] std::string toString(const std::string& var = "x") const {
        std::string sNum = numerator_.toString(var);
        std::string sDen = denominator_.toString(var);

        if (sNum.empty()) sNum = "0";

        // 分母が 1 なら分子のみ
        if (denominator_ == T(1)) {
            return sNum;
        }

        // 複数項なら括弧で囲む
        auto needsParens = [](const std::string& s) {
            return s.find('+') != std::string::npos ||
                   (s.find('-') != std::string::npos && s[0] != '-');
        };

        if (needsParens(sNum)) sNum = "(" + sNum + ")";
        if (needsParens(sDen)) sDen = "(" + sDen + ")";

        return sNum + "/" + sDen;
    }

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

    // ============================================================
    // 微小値切り捨て (浮動小数点型のみ)
    // ============================================================

    /// 絶対値が eps 未満の係数を 0 にし、その後約分を試みる
    void chopSmallValues(const T& eps) {
        static_assert(std::is_floating_point_v<T>,
                      "chopSmallValues requires floating-point coefficient type");
        auto chop = [&](Polynomial<T>& p) {
            bool changed = false;
            for (int i = 0; i <= p.degree(); ++i) {
                if (std::abs(p[i]) < eps) {
                    p[i] = T(0);
                    changed = true;
                }
            }
            if (changed) {
                // normalize: Polynomial(vector) 経由で再構築
                p = Polynomial<T>(std::vector<T>(
                    p.coefficients().begin(), p.coefficients().end()));
            }
        };
        chop(numerator_);
        chop(denominator_);
        simplify();
    }

    // ============================================================
    // 逆数
    // ============================================================

    /// Q(x)/P(x) を返す
    [[nodiscard]] RationalFunction reciprocal() const {
        return RationalFunction(denominator_, numerator_, false);
    }

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

    /// d/dx [P(x)/Q(x)] = (P'Q - PQ') / Q^2
    [[nodiscard]] RationalFunction derivative() const {
        return RationalFunction(
            numerator_.derivative() * denominator_
                - numerator_ * denominator_.derivative(),
            denominator_.pow(2)
        );
    }

private:
    // 分母の主係数が負なら、分子・分母の符号を反転して正規化
    void normalizeSign_() {
        if (denominator_.isZero()) return;  // 0 除算は呼び出し側の責任
        T lc = denominator_.leadingCoefficient();
        // 浮動小数点型の場合のみ符号を正規化
        if constexpr (std::is_floating_point_v<T>) {
            if (lc < T(0)) {
                numerator_ = -numerator_;
                denominator_ = -denominator_;
            }
        }
    }
};

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

template<PolynomialCoeff T>
[[nodiscard]] bool isZero(const RationalFunction<T>& f) {
    return f.isZero();
}

template<PolynomialCoeff T>
[[nodiscard]] bool isNonZero(const RationalFunction<T>& f) {
    return !f.isZero();
}

// ================================================================
// 合成: f(g(x)) — 有理関数 f に有理関数 g を代入
// ================================================================

/// f(g(x)) を計算する。
/// f = P(x)/Q(x), g = A(x)/B(x) のとき
/// f(g) = P(A/B) / Q(A/B)
///      = [Σ p_k A^k B^(K-k)] / [Σ q_k A^k B^(K-k)]
/// ただし K = max(deg P, deg Q)
template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> compose(
    const RationalFunction<T>& f,
    const RationalFunction<T>& g)
{
    const auto& P = f.numerator();
    const auto& Q = f.denominator();
    const auto& A = g.numerator();
    const auto& B = g.denominator();

    int K = std::max(P.degree(), Q.degree());
    if (K < 0) return RationalFunction<T>(T(0));

    // A^k と B^(K-k) を事前計算
    std::vector<Polynomial<T>> Apow(K + 1);
    std::vector<Polynomial<T>> Bpow(K + 1);
    Apow[0] = Polynomial<T>(T(1));
    Bpow[0] = Polynomial<T>(T(1));
    for (int i = 1; i <= K; ++i) {
        Apow[i] = Apow[i - 1] * A;
        Bpow[i] = Bpow[i - 1] * B;
    }

    // 分子: Σ p_k * A^k * B^(K-k)
    Polynomial<T> num;
    for (int k = 0; k <= P.degree(); ++k)
        num = num + Apow[k] * Bpow[K - k] * P[k];

    // 分母: Σ q_k * A^k * B^(K-k)
    Polynomial<T> den;
    for (int k = 0; k <= Q.degree(); ++k)
        den = den + Apow[k] * Bpow[K - k] * Q[k];

    return RationalFunction<T>(std::move(num), std::move(den));
}

// ================================================================
// 微小値切り捨て (自由関数版)
// ================================================================

template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> chopSmallValues(
    const RationalFunction<T>& f, const T& eps)
{
    RationalFunction<T> result = f;
    result.chopSmallValues(eps);
    return result;
}

// ================================================================
// 部分分数展開
// ================================================================
// partialFractions は polynomial_roots.hpp に依存するため、
// 使用する場合は RationalFunction_partialFractions.hpp を
// 別途 include すること。
// (循環依存を避けるためヘッダーを分離)

// ================================================================
// 型特性
// ================================================================
// 最小実現 (Minimal Realization)
// ================================================================

/// IIR フィルタ等の伝達関数 H(z) = B(z)/A(z) から共通因子を除去し、
/// 最小次数の等価な伝達関数を返す。
///
/// 用途: IIR フィルタ設計で生じる極零相殺 (pole-zero cancellation)
///
/// @param tf        伝達関数 (有理関数)
/// @param tolerance 近似 GCD の許容誤差
/// @return 最小実現された伝達関数
template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> minimalRealization(
    RationalFunction<T> tf,
    double tolerance = 1e-10) {

    tf.approximateSimplify(tolerance);
    return tf;
}

/// 分子・分母多項式から直接最小実現を行う。
///
/// @param numerator   分子多項式 B(z)
/// @param denominator 分母多項式 A(z)
/// @param tolerance   近似 GCD の許容誤差
/// @return 最小実現された伝達関数
template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> minimalRealization(
    const Polynomial<T>& numerator,
    const Polynomial<T>& denominator,
    double tolerance = 1e-10) {

    // doSimplify=false で構築し、近似約分を適用
    RationalFunction<T> tf(numerator, denominator, false);
    tf.approximateSimplify(tolerance);
    return tf;
}

// ================================================================

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

template<PolynomialCoeff T>
struct is_rational_function<RationalFunction<T>> : std::true_type {};

template<typename T>
inline constexpr bool is_rational_function_v = is_rational_function<T>::value;

} // namespace calx
