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

// RationalFunction.hpp
// Rational function template class (quotient of polynomials P(x)/Q(x))
//
// Design:
//   - Numerator and denominator are stored as Polynomial<T>
//   - Usable with T = double, float, Int, Float, Rational, Complex<T>, etc.
//   - Automatic cancellation via GCD
//   - O(n) evaluation using Horner's method

#pragma once

#include <math/core/Polynomial.hpp>

namespace sangi {

// ================================================================
// is_rational_function trait
// ================================================================
// NOTE: this must appear *before* the class definition. The requires-clause of
//   RationalFunction's member operator() uses is_rational_function_v *from within
//   the class*. Under the two-phase name lookup of g++/clang the name must be
//   declared before use (MSVC accepts it via lazy lookup of dependent names, so
//   the issue was previously hidden, but standalone g++ builds of the signal/IIR/
//   CAS headers that include this header indirectly all failed here). A partial
//   specialization only needs the type pattern, not a complete type, so a forward
//   declaration is sufficient.
template<PolynomialCoeff T> class RationalFunction;

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;

// ================================================================
// RationalFunction<T> class
// ================================================================

template<PolynomialCoeff T>
class RationalFunction {
private:
    Polynomial<T> numerator_;    // Numerator
    Polynomial<T> denominator_;  // Denominator

public:
    using value_type = T;

    // ============================================================
    // Constructors
    // ============================================================

    // Default: 0/1
    RationalFunction() : numerator_(T(0)), denominator_(T(1)) {}

    // From a scalar: c/1
    RationalFunction(const T& c) : numerator_(c), denominator_(T(1)) {}

    // From a polynomial: P(x)/1
    RationalFunction(const Polynomial<T>& p)
        : numerator_(p), denominator_(T(1)) {}

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

    // Specify numerator and denominator
    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();
    }

    // From initializer lists (ascending powers -- same convention as Polynomial)
    // Example: RationalFunction<double>({1,2,3}, {0,1}) -> (1+2x+3x^2)/(x)
    RationalFunction(std::initializer_list<T> num,
                     std::initializer_list<T> den,
                     bool doSimplify = true)
        : numerator_(num), denominator_(den)
    {
        normalizeSign_();
        if (doSimplify) simplify();
    }

    // ============================================================
    // Accessors
    // ============================================================

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

    // Numerator degree
    [[nodiscard]] int numeratorDegree() const { return numerator_.degree(); }

    // Denominator degree
    [[nodiscard]] int denominatorDegree() const { return denominator_.degree(); }

    // Whether the value is zero (numerator is zero)
    [[nodiscard]] bool isZero() const { return numerator_.isZero(); }

    // ============================================================
    // Cancellation (simplify numerator and denominator by 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_.quo(g);
            denominator_ = denominator_.quo(g);
            return true;
        }
        return false;
    }

    /// Approximate cancellation (simplify numerator and denominator using an
    /// approximate GCD computed by the epsilon-Euclidean method).
    ///
    /// Used when the floating-point coefficients of the rational function make
    /// the exact GCD collapse to 1 due to rounding error. Useful for minimal
    /// realization of IIR filters -- degree reduction via pole-zero cancellation.
    ///
    /// @param tolerance Maximum absolute value of remainder coefficients
    ///                  below which the remainder is treated as zero
    /// @return true if cancellation was performed
    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_.quo(g);
            denominator_ = denominator_.quo(g);
            // Renormalize because division rounding may leave tiny high-order terms
            normalizeSign_();
            return true;
        }
        return false;
    }

    // ============================================================
    // Evaluation (Horner's method)
    // ============================================================

    // Evaluate at a scalar value: compute P(x) / Q(x)
    template<typename U = T>
        requires (!is_polynomial_v<U> && !is_rational_function_v<U>)
    [[nodiscard]] U operator()(const U& x) const {
        return numerator_(x) / denominator_(x);
    }

    // Composition with a polynomial: P(q(x)) / Q(q(x)) -> returns a RationalFunction
    template<typename U>
        requires is_polynomial_v<U>
    [[nodiscard]] RationalFunction operator()(const U& x) const {
        return RationalFunction(numerator_(x), denominator_(x));
    }

    // Composition with a rational function: P(r(x)) / Q(r(x)) -> returns a RationalFunction
    [[nodiscard]] RationalFunction operator()(const RationalFunction& x) const {
        // numerator_(x), denominator_(x) return a RationalFunction via
        // Polynomial::operator() (Horner's method)
        RationalFunction num_result = numerator_(x);
        RationalFunction den_result = denominator_(x);
        return num_result / den_result;
    }

    // ============================================================
    // Comparison operators
    // ============================================================

    [[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);
    }

    // ============================================================
    // Unary operators
    // ============================================================

    [[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 + scalar
    // ============================================================

    [[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_);
    }

    // ============================================================
    // Compound assignment operators
    // ============================================================

    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; }

    // ============================================================
    // Exponentiation
    // ============================================================

    [[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: exponentiation of the reciprocal
        return RationalFunction(denominator_.pow(-n), numerator_.pow(-n), false);
    }

    // ============================================================
    // String conversion and stream output
    // ============================================================

    [[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";

        // If the denominator is 1, just emit the numerator
        if (denominator_ == T(1)) {
            return sNum;
        }

        // Wrap multi-term expressions in parentheses
        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();
    }

    // ============================================================
    // Chop small values (floating-point types only)
    // ============================================================

    /// Set coefficients whose absolute value is below eps to 0, then attempt to simplify
    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: reconstruct via Polynomial(vector)
                p = Polynomial<T>(std::vector<T>(
                    p.coefficients().begin(), p.coefficients().end()));
            }
        };
        chop(numerator_);
        chop(denominator_);
        simplify();
    }

    // ============================================================
    // Reciprocal
    // ============================================================

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

    // ============================================================
    // Differentiation
    // ============================================================

    /// 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:
    // Normalize signs: if the leading coefficient of the denominator is negative,
    // flip the signs of both numerator and denominator
    void normalizeSign_() {
        if (denominator_.isZero()) return;  // Division by zero is the caller's responsibility
        T lc = denominator_.leadingCoefficient();
        // Only normalize the sign for floating-point types
        if constexpr (std::is_floating_point_v<T>) {
            if (lc < T(0)) {
                numerator_ = -numerator_;
                denominator_ = -denominator_;
            }
        }
    }
};

// ================================================================
// Free functions
// ================================================================

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();
}

// ================================================================
// Composition: f(g(x)) -- substitute the rational function g into f
// ================================================================

/// Compute f(g(x)).
/// When f = P(x)/Q(x) and g = A(x)/B(x),
/// f(g) = P(A/B) / Q(A/B)
///      = [sum p_k A^k B^(K-k)] / [sum q_k A^k B^(K-k)]
/// where 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));

    // Precompute A^k and 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;
    }

    // Numerator: sum 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];

    // Denominator: sum 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));
}

// ================================================================
// Chop small values (free-function version)
// ================================================================

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

// ================================================================
// Partial fraction expansion
// ================================================================
// partialFractions depends on polynomial_roots.hpp, so to use it
// include RationalFunction_partialFractions.hpp separately.
// (The header is split to avoid a circular dependency.)

// ================================================================
// Type traits
// ================================================================
// Minimal realization
// ================================================================

/// Remove common factors from a transfer function H(z) = B(z)/A(z) (as used
/// for IIR filters) and return the equivalent minimum-degree transfer function.
///
/// Use case: pole-zero cancellation that arises in IIR filter design
///
/// @param tf        Transfer function (rational function)
/// @param tolerance Tolerance for the approximate GCD
/// @return Minimum-realized transfer function
template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> minimalRealization(
    RationalFunction<T> tf,
    double tolerance = 1e-10) {

    tf.approximateSimplify(tolerance);
    return tf;
}

/// Perform minimal realization directly from numerator and denominator polynomials.
///
/// @param numerator   Numerator polynomial B(z)
/// @param denominator Denominator polynomial A(z)
/// @param tolerance   Tolerance for the approximate GCD
/// @return Minimum-realized transfer function
template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> minimalRealization(
    const Polynomial<T>& numerator,
    const Polynomial<T>& denominator,
    double tolerance = 1e-10) {

    // Construct with doSimplify=false, then apply approximate cancellation
    RationalFunction<T> tf(numerator, denominator, false);
    tf.approximateSimplify(tolerance);
    return tf;
}

// ================================================================
// (the is_rational_function trait was moved above the class definition — for two-phase name lookup)
// ================================================================
// SANGI_POLY_DIV_RATIONAL mode
// ================================================================
//
// When SANGI_POLY_DIV_RATIONAL is defined, Polynomial / Polynomial
// returns a RationalFunction.
//
//   #define SANGI_POLY_DIV_RATIONAL
//   #include <math/core/RationalFunction.hpp>
//
//   Polynomial<double> p = {1, 0, 1};  // x^2 + 1
//   Polynomial<double> q = {1, 1};     // x + 1
//   auto rf = p / q;  // RationalFunction<double>: (x^2+1)/(x+1)
//
// Use p.quo(q) when you need the polynomial-division quotient instead.
// CMake: option(SANGI_POLY_DIV_RATIONAL ...) + add_compile_definitions(...)

#ifdef SANGI_POLY_DIV_RATIONAL
template<PolynomialCoeff T>
[[nodiscard]] RationalFunction<T> operator/(const Polynomial<T>& lhs, const Polynomial<T>& rhs) {
    return RationalFunction<T>(lhs, rhs);
}
#endif

} // namespace sangi