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

// approximation.hpp
// Function approximation algorithms
//
// Function list:
//   linearRegression(x, y)          → LinearRegressionResult
//   linearRegression(y)             → LinearRegressionResult (index-based)
//   weightedLinearRegression(x,y,w) → LinearRegressionResult
//   chebyshevApprox(f, a, b, n)     → ChebyshevApprox<T>  (Chebyshev approximation)
//   padeApprox(taylor, M, N)        → PadeResult<T>       (Padé approximation)
//   padeTable(taylor, Mmax, Nmax)   → vector<vector<PadeResult<T>>>       (Padé table)
//   taylorToSFraction(taylor, kmax) → StieltjesFraction<T>  (qd-algorithm)
//   taylorToJFraction(taylor, kmax) → JacobiFraction<T>      (via qd-algorithm)
//   evalSFraction(sf, x), evalJFraction(jf, x)               (continued-fraction evaluation)
//   rationalInterpolation(x, y, M, N) → RationalInterpResult<T> (multipoint Padé / Cauchy)
//   evalRational(result, x)                                  (rational function evaluation)
//   matrixExpPade(A, m=6)           → Matrix<T>           (matrix exponential e^A, Scaling-and-Squaring)
//   polynomialFit(x, y, degree)     → std::vector<T>      (least-squares polynomial approximation)
//
// Selection guide:
//   - Priority on compute cost / easy coefficient analysis: chebyshevApprox
//   - Taylor coefficients known + accuracy improvement: padeApprox
//   - Computing e^A: matrixExpPade (Scaling-and-Squaring, m=6 is standard accuracy)
//   - Fitting a polynomial to data points: polynomialFit (note: high degree worsens the condition number)

#ifndef SANGI_APPROXIMATION_HPP
#define SANGI_APPROXIMATION_HPP

#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <math/core/generic_math.hpp>
#include <math/concepts/algebraic_concepts.hpp>
#include <math/linalg/solvers.hpp>
#include <math/linalg/decomposition.hpp>
#include <math/linalg/matrix_functions.hpp>   // algorithms::expm (matrixExpPade delegates here)
#include <cmath>
#include <functional>
#include <span>
#include <stdexcept>
#include <vector>

namespace sangi {

// =====================================================================
// Result structures
// =====================================================================

/**
 * @brief Result of linear regression
 * @note slope/intercept/R² are now made explicit
 */
template<concepts::OrderedField T>
struct LinearRegressionResult {
    T slope;          // slope a
    T intercept;      // intercept b (y = a*x + b)
    T r_squared;      // coefficient of determination R²
};

// =====================================================================
// Linear regression
// =====================================================================

/**
 * @brief Linear regression y = a*x + b
 * @param x array of independent variables
 * @param y array of dependent variables (same size as x)
 * @return LinearRegressionResult{slope, intercept, r_squared}
 * @throw std::invalid_argument on size mismatch, or fewer than 2 elements
 */
template<concepts::OrderedField T>
[[nodiscard]] LinearRegressionResult<T> linearRegression(
    std::span<const T> x, std::span<const T> y) {

    if (x.size() != y.size()) {
        throw std::invalid_argument("linearRegression: x and y must have the same size");
    }
    if (x.size() < 2) {
        throw std::invalid_argument("linearRegression: need at least 2 data points");
    }

    std::size_t n = x.size();
    T n_val = static_cast<T>(n);

    // means
    T mx = T(0), my = T(0);
    for (std::size_t i = 0; i < n; ++i) {
        mx += x[i];
        my += y[i];
    }
    mx /= n_val;
    my /= n_val;

    // Sxx, Sxy, Syy
    T sxx = T(0), sxy = T(0), syy = T(0);
    for (std::size_t i = 0; i < n; ++i) {
        T dx = x[i] - mx;
        T dy = y[i] - my;
        sxx += dx * dx;
        sxy += dx * dy;
        syy += dy * dy;
    }

    if (sxx == T(0)) {
        throw std::invalid_argument("linearRegression: all x values are identical");
    }

    T slope = sxy / sxx;
    T intercept = my - slope * mx;

    // R² = Sxy² / (Sxx * Syy)
    T r_squared = (syy == T(0)) ? T(1) : (sxy * sxy) / (sxx * syy);

    return {slope, intercept, r_squared};
}

template<concepts::OrderedField T>
[[nodiscard]] LinearRegressionResult<T> linearRegression(
    const std::vector<T>& x, const std::vector<T>& y) {
    return linearRegression(std::span<const T>(x), std::span<const T>(y));
}

/**
 * @brief Index-based linear regression y = a*i + b (i = 0, 1, ..., n-1)
 * @param y array of dependent variables
 * @return LinearRegressionResult{slope, intercept, r_squared}
 */
template<concepts::OrderedField T>
[[nodiscard]] LinearRegressionResult<T> linearRegression(std::span<const T> y) {
    if (y.size() < 2) {
        throw std::invalid_argument("linearRegression: need at least 2 data points");
    }
    std::vector<T> x(y.size());
    for (std::size_t i = 0; i < y.size(); ++i) {
        x[i] = static_cast<T>(i);
    }
    return linearRegression(std::span<const T>(x), y);
}

template<concepts::OrderedField T>
[[nodiscard]] LinearRegressionResult<T> linearRegression(const std::vector<T>& y) {
    return linearRegression(std::span<const T>(y));
}

/**
 * @brief Weighted linear regression y = a*x + b
 * @param x independent variable
 * @param y dependent variable
 * @param weights weights
 * @return LinearRegressionResult{slope, intercept, r_squared}
 * @note stabilized with algorithms::solve
 */
template<concepts::OrderedField T>
[[nodiscard]] LinearRegressionResult<T> weightedLinearRegression(
    std::span<const T> x, std::span<const T> y, std::span<const T> weights) {

    if (x.size() != y.size() || x.size() != weights.size()) {
        throw std::invalid_argument("weightedLinearRegression: x, y, weights must have the same size");
    }
    if (x.size() < 2) {
        throw std::invalid_argument("weightedLinearRegression: need at least 2 data points");
    }

    // Normal equations: [Σw*x², Σw*x; Σw*x, Σw] [a; b] = [Σw*x*y; Σw*y]
    T wxx = T(0), wx = T(0), w_sum = T(0), wxy = T(0), wy = T(0);
    for (std::size_t i = 0; i < x.size(); ++i) {
        T wi = weights[i];
        wxx += wi * x[i] * x[i];
        wx  += wi * x[i];
        w_sum += wi;
        wxy += wi * x[i] * y[i];
        wy  += wi * y[i];
    }

    if (w_sum == T(0)) {
        throw std::invalid_argument("weightedLinearRegression: sum of weights must not be zero");
    }

    // Solve the 2x2 linear system with algorithms::solve
    Matrix<T> A(2, 2);
    A(0, 0) = wxx;  A(0, 1) = wx;
    A(1, 0) = wx;   A(1, 1) = w_sum;

    Vector<T> b(2);
    b[0] = wxy;
    b[1] = wy;

    Vector<T> solution = algorithms::solve(A, b, algorithms::SolverType::LU);

    T slope = solution[0];
    T intercept = solution[1];

    // weighted R²
    T wmx = wx / w_sum;
    T wmy = wy / w_sum;
    T ss_tot = T(0), ss_res = T(0);
    for (std::size_t i = 0; i < x.size(); ++i) {
        T wi = weights[i];
        T y_pred = slope * x[i] + intercept;
        ss_res += wi * (y[i] - y_pred) * (y[i] - y_pred);
        ss_tot += wi * (y[i] - wmy) * (y[i] - wmy);
    }
    T r_squared = (ss_tot == T(0)) ? T(1) : T(1) - ss_res / ss_tot;

    return {slope, intercept, r_squared};
}

template<concepts::OrderedField T>
[[nodiscard]] LinearRegressionResult<T> weightedLinearRegression(
    const std::vector<T>& x, const std::vector<T>& y, const std::vector<T>& weights) {
    return weightedLinearRegression(
        std::span<const T>(x), std::span<const T>(y), std::span<const T>(weights));
}

// =====================================================================
// Chebyshev approximation
// =====================================================================

/**
 * @brief Class holding the result of a Chebyshev approximation
 *
 * Approximates a function f by a degree-n Chebyshev polynomial on [a, b].
 * The Chebyshev coefficients are computed at construction, and operator() evaluates the approximation.
 *
 * @note Improvements:
 *   - class template now has concepts constraints
 *   - function pointer → std::function
 *   - assert → exception
 *   - toCoefficients() can convert to ordinary polynomial coefficients
 */
template<concepts::OrderedField T>
class ChebyshevApprox {
public:
    /**
     * @brief Construct a Chebyshev approximation
     * @param f function to approximate
     * @param a left endpoint of the interval
     * @param b right endpoint of the interval
     * @param n degree of the Chebyshev polynomial (≥ 0)
     */
    ChebyshevApprox(const std::function<T(T)>& f, T a, T b, std::size_t n)
        : m_a(a), m_b(b), m_n(n + 1), m_c(n + 1) {

        if (a >= b) {
            throw std::invalid_argument("ChebyshevApprox: a must be less than b");
        }

        const T pi = std::acos(T(-1));
        const T bma = T(0.5) * (b - a);
        const T bpa = T(0.5) * (b + a);
        const std::size_t N = m_n;

        // Function values at the Chebyshev nodes
        std::vector<T> fvals(N);
        for (std::size_t k = 0; k < N; ++k) {
            T y = std::cos(pi * (static_cast<T>(k) + T(0.5)) / static_cast<T>(N));
            fvals[k] = f(y * bma + bpa);
        }

        // Compute the Chebyshev coefficients (DCT-II)
        T fac = T(2) / static_cast<T>(N);
        for (std::size_t j = 0; j < N; ++j) {
            T sum = T(0);
            for (std::size_t k = 0; k < N; ++k) {
                sum += fvals[k] * std::cos(pi * static_cast<T>(j) * (static_cast<T>(k) + T(0.5)) / static_cast<T>(N));
            }
            m_c[j] = fac * sum;
        }
    }

    /**
     * @brief Evaluate the approximation (Clenshaw recurrence)
     * @param x evaluation point (inside the interval [a, b])
     * @return approximate value
     */
    [[nodiscard]] T operator()(T x) const {
        return evaluate(x, m_n);
    }

    /**
     * @brief Evaluate the approximation using only the first m terms (truncation)
     * @param x evaluation point
     * @param m number of terms to use (≤ m_n)
     * @return approximate value
     */
    [[nodiscard]] T evaluate(T x, std::size_t m) const {
        if (m > m_n) m = m_n;
        if (m == 0) return T(0);

        // Map x to [-1, 1]
        T y = (T(2) * x - m_a - m_b) / (m_b - m_a);

        // Clenshaw recurrence
        T d = T(0), dd = T(0);
        T y2 = T(2) * y;
        for (std::size_t j = m - 1; j >= 1; --j) {
            T sv = d;
            d = y2 * d - dd + m_c[j];
            dd = sv;
        }
        return y * d - dd + T(0.5) * m_c[0];
    }

    /** @brief Get the Chebyshev coefficients */
    [[nodiscard]] const std::vector<T>& coefficients() const { return m_c; }

    /** @brief Get the number of terms */
    [[nodiscard]] std::size_t size() const { return m_n; }

    /** @brief Get the interval */
    [[nodiscard]] T lower() const { return m_a; }
    [[nodiscard]] T upper() const { return m_b; }

    /**
     * @brief Convert to ordinary polynomial coefficients [a0, a1, ..., aN]
     *
     * Computes the coefficients of p(x) = Σ a_j * x^j from the Chebyshev coefficients c_k.
     * Includes the change of variable [a,b] → [-1,1].
     * @return polynomial coefficients (a0 + a1*x + ... + aN*x^N)
     */
    [[nodiscard]] std::vector<T> toPolynomialCoefficients() const {
        // Convert the Chebyshev polynomial expansion to an ordinary polynomial
        // T_0(y) = 1, T_1(y) = y, T_{k+1}(y) = 2*y*T_k(y) - T_{k-1}(y)
        // where y = (2*x - a - b) / (b - a)

        const std::size_t N = m_n;
        // First build the polynomial coefficients on [-1,1]
        // d[k][j] = degree-j coefficient of T_k(y)
        std::vector<std::vector<T>> tcheb(N);
        tcheb[0].assign(1, T(1));  // T_0 = 1
        if (N > 1) {
            tcheb[1].assign(2, T(0));
            tcheb[1][1] = T(1);  // T_1 = y
        }
        for (std::size_t k = 2; k < N; ++k) {
            tcheb[k].assign(k + 1, T(0));
            // T_k = 2*y*T_{k-1} - T_{k-2}
            for (std::size_t j = 0; j < tcheb[k - 1].size(); ++j) {
                if (j + 1 < tcheb[k].size())
                    tcheb[k][j + 1] += T(2) * tcheb[k - 1][j];
            }
            for (std::size_t j = 0; j < tcheb[k - 2].size(); ++j) {
                tcheb[k][j] -= tcheb[k - 2][j];
            }
        }

        // Weighted sum with the Chebyshev coefficients (c[0]/2 * T_0 + c[1]*T_1 + ...)
        std::vector<T> poly_y(N, T(0));
        for (std::size_t k = 0; k < N; ++k) {
            T weight = (k == 0) ? m_c[0] / T(2) : m_c[k];
            for (std::size_t j = 0; j < tcheb[k].size(); ++j) {
                poly_y[j] += weight * tcheb[k][j];
            }
        }

        // Apply the change of variable y = (2*x - a - b) / (b - a)
        // y = alpha * x + beta, alpha = 2/(b-a), beta = -(a+b)/(b-a)
        T alpha = T(2) / (m_b - m_a);
        T beta = -(m_a + m_b) / (m_b - m_a);

        // Substitute (alpha*x + beta) for y in p(y)
        // Result: poly_x[j] = coefficient of x^j
        std::vector<T> poly_x(N, T(0));
        // Compute the coefficients of (alpha*x + beta)^k by binomial expansion
        // powers_of_ax_plus_b[k][j] = degree-j coefficient of (alpha*x + beta)^k
        std::vector<T> current(1, T(1));  // (alpha*x + beta)^0 = 1
        for (std::size_t k = 0; k < N; ++k) {
            // Add poly_y[k] * current into poly_x
            for (std::size_t j = 0; j < current.size(); ++j) {
                poly_x[j] += poly_y[k] * current[j];
            }
            // current *= (alpha*x + beta)
            std::vector<T> next(current.size() + 1, T(0));
            for (std::size_t j = 0; j < current.size(); ++j) {
                next[j] += beta * current[j];
                next[j + 1] += alpha * current[j];
            }
            current = std::move(next);
        }

        return poly_x;
    }

private:
    T m_a, m_b;
    std::size_t m_n;    // number of terms (degree + 1)
    std::vector<T> m_c; // Chebyshev coefficients
};

// =====================================================================
// Padé approximation
// =====================================================================

/**
 * @brief Result of a Padé approximation
 * @note Form p(x)/q(x). q[0] = 1 (normalized)
 */
template<concepts::OrderedField T>
struct PadeResult {
    std::vector<T> numerator;    // numerator polynomial coefficients [p0, p1, ..., pM]
    std::vector<T> denominator;  // denominator polynomial coefficients [1, q1, ..., qN]
    bool valid;                  // whether a solution was obtained
};

/**
 * @brief Padé approximation — build a rational-function approximation from a Taylor series
 *
 * Computes the [M/N] Padé approximation from Taylor coefficients a[0], a[1], ..., a[M+N].
 * The result is p(x)/q(x), normalized so that q(0) = 1.
 *
 * @param taylor Taylor expansion coefficients (at least M+N+1 needed)
 * @param M degree of the numerator polynomial
 * @param N degree of the denominator polynomial
 * @return PadeResult{numerator, denominator, valid}
 *
 * @note Improvements:
 *   - removed cout debug output
 *   - stable linear-system solving with algorithms::solve
 *   - self-contained, no dependency on a rational-expression class
 */
template<concepts::OrderedField T>
[[nodiscard]] PadeResult<T> padeApprox(
    std::span<const T> taylor,
    std::size_t M, std::size_t N) {

    const std::size_t L = M + N;
    if (taylor.size() < L + 1) {
        throw std::invalid_argument("padeApprox: need at least M+N+1 Taylor coefficients");
    }

    // Linear system: L+1 unknowns [p0, ..., pM, q1, ..., qN]
    // Condition: Σ_{k=0}^{M} p_k * x^k = (Σ_{j=0}^{L} a_j * x^j) * (1 + Σ_{l=1}^{N} q_l * x^l) mod x^{L+1}
    // Expanding this, at each degree m:
    //   m ≤ M: p_m + Σ_{l=1}^{min(m,N)} q_l * a_{m-l} = a_m
    //   m > M: Σ_{l=1}^{min(m,N)} q_l * a_{m-l} = a_m  (where m-l ≥ 0)

    // First solve for q (a linear system in N unknowns)
    if (N == 0) {
        // denominator = 1 → numerator = the Taylor coefficients as-is
        std::vector<T> num(taylor.begin(), taylor.begin() + static_cast<std::ptrdiff_t>(M + 1));
        return {num, {T(1)}, true};
    }

    // Linear system for q: for degrees M+1, M+2, ..., M+N
    // a_{m} + q_1 * a_{m-1} + q_2 * a_{m-2} + ... + q_N * a_{m-N} = 0  (m = M+1..M+N)
    Matrix<T> A(static_cast<typename Matrix<T>::size_type>(N),
                 static_cast<typename Matrix<T>::size_type>(N));
    Vector<T> b(static_cast<typename Vector<T>::size_type>(N));

    for (std::size_t i = 0; i < N; ++i) {
        std::size_t m = M + 1 + i;
        for (std::size_t j = 0; j < N; ++j) {
            std::size_t idx = m - 1 - j;  // a_{m-(j+1)}
            A(static_cast<typename Matrix<T>::size_type>(i),
              static_cast<typename Matrix<T>::size_type>(j)) =
                (idx <= L) ? taylor[idx] : T(0);
        }
        b[static_cast<typename Vector<T>::size_type>(i)] = -taylor[m];
    }

    // Solve for q
    Vector<T> q_vec(static_cast<typename Vector<T>::size_type>(0));
    try {
        q_vec = algorithms::solve(A, b, algorithms::SolverType::LU);
    } catch (...) {
        return {{}, {}, false};
    }

    // Denominator coefficients: [1, q1, q2, ..., qN]
    std::vector<T> denom(N + 1);
    denom[0] = T(1);
    for (std::size_t j = 0; j < N; ++j) {
        denom[j + 1] = q_vec[static_cast<typename Vector<T>::size_type>(j)];
    }

    // Numerator coefficients: p_m = a_m + q_1*a_{m-1} + ... + q_min(m,N)*a_{m-min(m,N)}
    std::vector<T> num(M + 1);
    for (std::size_t m = 0; m <= M; ++m) {
        T val = taylor[m];
        for (std::size_t l = 1; l <= std::min(m, N); ++l) {
            val += denom[l] * taylor[m - l];
        }
        num[m] = val;
    }

    return {num, denom, true};
}

template<concepts::OrderedField T>
[[nodiscard]] PadeResult<T> padeApprox(
    const std::vector<T>& taylor, std::size_t M, std::size_t N) {
    return padeApprox(std::span<const T>(taylor), M, N);
}

/**
 * @brief Evaluate a Padé approximation
 * @param result result of padeApprox
 * @param x evaluation point
 * @return p(x) / q(x)
 */
template<concepts::OrderedField T>
[[nodiscard]] T evaluatePade(const PadeResult<T>& result, T x) {
    // Evaluate numerator and denominator by Horner's method
    T num = result.numerator.back();
    for (std::size_t i = result.numerator.size() - 1; i > 0; --i) {
        num = num * x + result.numerator[i - 1];
    }

    T den = result.denominator.back();
    for (std::size_t i = result.denominator.size() - 1; i > 0; --i) {
        den = den * x + result.denominator[i - 1];
    }

    return num / den;
}

// =====================================================================
// Padé table — generate the [M/N] triangular array in one pass
// =====================================================================

/**
 * @brief Padé table — compute the Padé approximation for all (M, N) with 0 ≤ M ≤ Mmax, 0 ≤ N ≤ Nmax
 *
 * From the input Taylor coefficients, returns [M/N] Padé as a table with rows M and columns N.
 * table[M][N] is a PadeResult<T>, laid out so that the diagonal [M/M] and antidiagonal [M+N=const]
 * are easy to compare.
 *
 * Degenerate cells (valid=false from ill-conditioning) are returned as-is with valid=false.
 *
 * @param taylor Taylor coefficients (at least Mmax + Nmax + 1)
 * @param Mmax   maximum numerator degree
 * @param Nmax   maximum denominator degree
 * @return a (Mmax+1) × (Nmax+1) two-dimensional array, table[m][n] = [m/n] Padé
 *
 * @note Complexity is O(Mmax · Nmax · N²) (each cell solves an N×N system).
 *       For large tables, when you only want e.g. the diagonal [M/M], it is often cheaper
 *       to call padeApprox directly.
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::vector<PadeResult<T>>>
padeTable(std::span<const T> taylor, std::size_t Mmax, std::size_t Nmax) {
    if (taylor.size() < Mmax + Nmax + 1)
        throw std::invalid_argument("padeTable: need at least Mmax+Nmax+1 Taylor coefficients");

    std::vector<std::vector<PadeResult<T>>> table(Mmax + 1);
    for (std::size_t m = 0; m <= Mmax; ++m) {
        table[m].reserve(Nmax + 1);
        for (std::size_t n = 0; n <= Nmax; ++n) {
            try {
                table[m].push_back(padeApprox<T>(taylor, m, n));
            } catch (...) {
                table[m].push_back({{}, {}, false});
            }
        }
    }
    return table;
}

template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::vector<PadeResult<T>>>
padeTable(const std::vector<T>& taylor, std::size_t Mmax, std::size_t Nmax) {
    return padeTable(std::span<const T>(taylor), Mmax, Nmax);
}

// =====================================================================
// Continued-fraction conversions — Stieltjes (S-fraction) / Jacobi (J-fraction)
// =====================================================================

/**
 * @brief Stieltjes continued fraction (S-fraction)
 *
 *   f(z) = c_0 / (1 + a_1 z / (1 + a_2 z / (1 + a_3 z / ...)))
 *
 * When terminated (terminated=true), a 0 appears somewhere in a and everything after is truncated.
 * A pure polynomial f is represented with a empty (only b0).
 */
template<concepts::OrderedField T>
struct StieltjesFraction {
    T b0;                  // c_0 (leading constant, usually taylor[0])
    std::vector<T> a;      // a_1, a_2, a_3, ... (sign-flipped qd coefficients)
    bool terminated;       // whether the qd table hit 0 partway and the continued fraction completed at finite length
    int order;             // number of input Taylor coefficients - 1
};

/**
 * @brief Jacobi continued fraction (J-fraction)
 *
 *   f(z) = c_0 / (1 - β_0 z - α_1² z² / (1 - β_1 z - α_2² z² / ...))
 *
 * Corresponds to the three-term recurrence of orthogonal polynomials. Convergence of the J-fraction
 * up to stage n equals the diagonal [n/n] Padé approximation. It groups two S-fraction stages at a time.
 */
template<concepts::OrderedField T>
struct JacobiFraction {
    T c0;                       // leading constant
    std::vector<T> beta;        // β_0, β_1, β_2, ...
    std::vector<T> alpha2;      // α_1², α_2², ... (one fewer than β, or the same count)
    bool terminated;
    int order;
};

namespace detail_qd {

/**
 * @brief Build the qd table from Taylor coefficients via the Rutishauser qd-algorithm
 *
 * From input c_0, c_1, ..., c_{2K} (= 2K+1 values), compute up to depth K:
 *   q_k^(0) for k=1..K
 *   e_k^(0) for k=1..K-1
 *
 * If degeneracy (q_k or e_k being 0) is detected, terminate at that point with termination=true.
 *
 * Return form: a pair (q_levels, e_levels), both vectors of length K, but elements after
 * termination are undefined. The actual length is conveyed via last_valid.
 */
template<concepts::OrderedField T>
inline void runQdAlgorithm(
    std::span<const T> c, int K,
    std::vector<T>& q0,    // output: q_1^(0), q_2^(0), ..., q_K^(0)
    std::vector<T>& e0,    // output: e_1^(0), e_2^(0), ..., e_{K-1}^(0)
    bool& terminated)
{
    terminated = false;
    q0.clear();
    e0.clear();

    if (K < 1) return;
    if (c[0] == T(0)) {
        // If the leading term is 0, q_1 = c_1 / c_0 cannot be defined
        terminated = true;
        return;
    }

    // For row n, q_1^(n) = c_{n+1}/c_n, e_0^(n) = 0.
    // We need q_1^(n) for n = 0..(2K-1); e_0 is not needed (always 0).
    const int Nrows = 2 * K;
    if (static_cast<int>(c.size()) < Nrows + 1) {
        // Shrink to the largest K computable with the given coefficients
        K = static_cast<int>(c.size() - 1) / 2;
        if (K < 1) { terminated = true; return; }
    }

    // q[k][n] = q_{k+1}^(n) (1-indexed → 0-indexed shift)
    // Dynamic programming: for column k = 1..K, keep rows n = 0..(2K - 2k + 1)
    //           q_1^(n) is n = 0..2K-1
    //           q_2^(n) is n = 0..2K-3
    //           ...
    //           q_K^(n) is n = 0..1
    //           e_k^(n) is n = 0..(2K - 2k)
    // For simplicity, instead of managing two rows at a time, use a 2D vector:
    std::vector<std::vector<T>> q(K), e(K);
    // q[k-1] = the q_k^(n) column, e[k-1] = the e_k^(n) column

    // k = 1: q_1^(n) = c_{n+1} / c_n
    int rows1 = std::min(2 * K, static_cast<int>(c.size()) - 1);
    q[0].resize(rows1);
    for (int n = 0; n < rows1; ++n) {
        if (c[n] == T(0)) { terminated = true; q[0].resize(n); break; }
        q[0][n] = c[n + 1] / c[n];
    }
    if (q[0].empty()) { terminated = true; return; }
    q0.push_back(q[0][0]);

    // Iterate: k = 2..K
    for (int k = 2; k <= K; ++k) {
        // e_{k-1}^(n) = q_{k-1}^(n+1) - q_{k-1}^(n) + e_{k-2}^(n+1)
        // (when k-2 = 0, e_0^(n) = 0)
        const int prev_q_size = static_cast<int>(q[k - 2].size());
        const int e_size = prev_q_size - 1;
        if (e_size < 1) { terminated = true; break; }
        e[k - 2].resize(e_size);
        for (int n = 0; n < e_size; ++n) {
            T prev_e = (k - 2 >= 1 && n + 1 < static_cast<int>(e[k - 3].size()))
                       ? e[k - 3][n + 1] : T(0);
            e[k - 2][n] = q[k - 2][n + 1] - q[k - 2][n] + prev_e;
        }
        if (k - 2 == 0 || k == 2) {
            // output e_1^(0)
        }
        e0.push_back(e[k - 2][0]);

        // Termination test: e_{k-1}^(0) == 0
        if (e[k - 2][0] == T(0)) { terminated = true; break; }

        // q_k^(n) = q_{k-1}^(n+1) · e_{k-1}^(n+1) / e_{k-1}^(n)
        const int q_size = e_size - 1;
        if (q_size < 1) { terminated = true; break; }
        q[k - 1].resize(q_size);
        for (int n = 0; n < q_size; ++n) {
            if (e[k - 2][n] == T(0)) { terminated = true; q[k - 1].resize(n); break; }
            q[k - 1][n] = q[k - 2][n + 1] * e[k - 2][n + 1] / e[k - 2][n];
        }
        if (q[k - 1].empty()) { terminated = true; break; }
        q0.push_back(q[k - 1][0]);

        // q_k^(0) == 0 also terminates (Padé degeneracy)
        if (q[k - 1][0] == T(0)) { terminated = true; break; }
    }
}

} // namespace detail_qd

/**
 * @brief Convert a Taylor series to a Stieltjes continued fraction
 *
 *   f(z) = c_0 + c_1 z + c_2 z² + ...
 *        = c_0 / (1 + a_1 z / (1 + a_2 z / ...))
 *
 * @param taylor Taylor coefficients c_0, c_1, ..., c_L
 * @param kmax   maximum number of a to generate (default: automatic based on input)
 *
 * @return StieltjesFraction { b0=c_0, a=[a_1, a_2, ...], terminated, order }
 *
 * Truncate when degeneracy (e_k or q_k being 0) is detected, terminated=true.
 *
 * Correspondence between S-fraction coefficients and the qd table:
 *   a_1 = -q_1^(0),  a_2 = -e_1^(0),  a_3 = -q_2^(0),  a_4 = -e_2^(0),  ...
 */
template<concepts::OrderedField T>
[[nodiscard]] StieltjesFraction<T>
taylorToSFraction(std::span<const T> taylor, int kmax = -1) {
    if (taylor.empty())
        throw std::invalid_argument("taylorToSFraction: empty taylor");

    const int L = static_cast<int>(taylor.size()) - 1;
    if (kmax < 0) kmax = L;
    // Coefficients needed by qd: 2K gives q_K^(0), 2K-1 gives e_{K-1}^(0)
    int K = std::min(kmax / 2 + 1, (L + 1) / 2);  // round conservatively
    if (K < 1) K = 1;

    StieltjesFraction<T> sf;
    sf.b0 = taylor[0];
    sf.terminated = false;
    sf.order = L;

    if (taylor[0] == T(0)) {
        // Leading 0 → b0=0, empty continued fraction (a shift etc. is expected at the caller level)
        sf.terminated = true;
        return sf;
    }

    std::vector<T> q0, e0;
    detail_qd::runQdAlgorithm<T>(taylor, K, q0, e0, sf.terminated);

    // a_{2k-1} = -q_k^(0), a_{2k} = -e_k^(0)
    sf.a.reserve(q0.size() + e0.size());
    for (std::size_t i = 0; i < std::max(q0.size(), e0.size() + 1); ++i) {
        if (i < q0.size()) {
            sf.a.push_back(-q0[i]);
            if (static_cast<int>(sf.a.size()) >= kmax) break;
        }
        if (i < e0.size()) {
            sf.a.push_back(-e0[i]);
            if (static_cast<int>(sf.a.size()) >= kmax) break;
        }
    }
    return sf;
}

template<concepts::OrderedField T>
[[nodiscard]] StieltjesFraction<T>
taylorToSFraction(const std::vector<T>& taylor, int kmax = -1) {
    return taylorToSFraction<T>(std::span<const T>(taylor), kmax);
}

/**
 * @brief Convert a Taylor series to a Jacobi continued fraction
 *
 *   f(z) = c_0 / (1 - β_0 z - α_1² z² / (1 - β_1 z - α_2² z² / ...))
 *
 * Groups two S-fraction stages at a time. Equivalent to the three-term recurrence of orthogonal polynomials.
 *
 * Conversion formulas (from the qd table):
 *   β_0 = q_1^(0)
 *   β_n = q_{n+1}^(0) + e_n^(0)   (n ≥ 1)
 *   α_n² = q_n^(0) · e_n^(0)        (n ≥ 1)
 */
template<concepts::OrderedField T>
[[nodiscard]] JacobiFraction<T>
taylorToJFraction(std::span<const T> taylor, int kmax = -1) {
    if (taylor.empty())
        throw std::invalid_argument("taylorToJFraction: empty taylor");

    const int L = static_cast<int>(taylor.size()) - 1;
    if (kmax < 0) kmax = (L + 1) / 2;
    int K = std::min(kmax + 1, (L + 1) / 2);
    if (K < 1) K = 1;

    JacobiFraction<T> jf;
    jf.c0 = taylor[0];
    jf.terminated = false;
    jf.order = L;

    if (taylor[0] == T(0)) { jf.terminated = true; return jf; }

    std::vector<T> q0, e0;
    detail_qd::runQdAlgorithm<T>(taylor, K, q0, e0, jf.terminated);

    // β_0 = q_1^(0)
    if (q0.empty()) return jf;
    jf.beta.push_back(q0[0]);

    // β_n = q_{n+1}^(0) + e_n^(0), α_n² = q_n^(0) · e_n^(0)
    for (std::size_t n = 0; n < e0.size(); ++n) {
        // α_{n+1}² = q_{n+1}^(0) · e_{n+1}^(0)?  Unify the convention: with α_n² (n ≥ 1),
        // the definition is α_n² = q_n^(0) · e_n^(0), i.e. it starts at α_1².
        T alpha_sq = q0[n] * e0[n];
        jf.alpha2.push_back(alpha_sq);
        if (n + 1 < q0.size()) {
            jf.beta.push_back(q0[n + 1] + e0[n]);
        }
        if (static_cast<int>(jf.beta.size()) > kmax) break;
    }
    return jf;
}

template<concepts::OrderedField T>
[[nodiscard]] JacobiFraction<T>
taylorToJFraction(const std::vector<T>& taylor, int kmax = -1) {
    return taylorToJFraction<T>(std::span<const T>(taylor), kmax);
}

/**
 * @brief Evaluate a Stieltjes continued fraction (bottom-up)
 *
 *   f(z) = b0 / (1 + a_1 z / (1 + a_2 z / (1 + a_K z)))
 */
template<concepts::OrderedField T>
[[nodiscard]] T evalSFraction(const StieltjesFraction<T>& sf, T x) {
    if (sf.a.empty()) return sf.b0;
    T acc = T(1) + sf.a.back() * x;
    for (int i = static_cast<int>(sf.a.size()) - 2; i >= 0; --i) {
        acc = T(1) + sf.a[i] * x / acc;
    }
    return sf.b0 / acc;
}

/**
 * @brief Evaluate a Jacobi continued fraction (bottom-up)
 *
 *   f(z) = c0 / (1 - β_0 z - α_1² z² / (1 - β_1 z - α_2² z² / (1 - β_K z)))
 *
 * The α_n² array is one shorter than β (the standard J-fraction form). If the lengths are equal,
 * use up to the trailing α² (accounting for truncation).
 */
template<concepts::OrderedField T>
[[nodiscard]] T evalJFraction(const JacobiFraction<T>& jf, T x) {
    if (jf.beta.empty()) return jf.c0;
    const int K = static_cast<int>(jf.beta.size());
    T acc = T(1) - jf.beta[K - 1] * x;
    for (int i = K - 2; i >= 0; --i) {
        T num = (i + 1 <= static_cast<int>(jf.alpha2.size())) ? jf.alpha2[i] : T(0);
        acc = T(1) - jf.beta[i] * x - num * x * x / acc;
    }
    return jf.c0 / acc;
}

// =====================================================================
// Multipoint Padé / rational interpolation (Cauchy interpolation)
// =====================================================================

/**
 * @brief Result of rational interpolation
 *
 * A generalization of standard Padé (fitting Taylor coefficients at one point) to multipoint
 * (fitting values at distinct x_i). Also called Cauchy / Newton-Padé / multipoint Padé.
 *
 * Returns the [M/N] rational function p(x)/q(x) satisfying f(x_i) = y_i (M + N + 1 = number of data points).
 * When q(x_i) = 0 (a pole), unless p(x_i) is simultaneously 0 (zero-pole cancellation), valid=false.
 */
template<concepts::OrderedField T>
struct RationalInterpResult {
    std::vector<T> numerator;     // [p0, p1, ..., pM]
    std::vector<T> denominator;   // [1, q1, ..., qN] (normalized to q0 = 1)
    bool valid;
};

/**
 * @brief Multipoint Padé / rational interpolation — build a rational function from values y_i at distinct x_i
 *
 * Linear system: for each i (i = 0..M+N)
 *   p(x_i) - y_i · q(x_i) = 0
 *   ⟺ Σ_{k=0..M} p_k x_i^k - Σ_{k=1..N} y_i q_k x_i^k = y_i
 * Solve the (M+N+1)-variable system.
 *
 * @param x data points (M+N+1, distinct)
 * @param y function values (same size as x)
 * @param M numerator degree
 * @param N denominator degree
 *
 * Note: if M + N + 1 ≠ x.size(), throws invalid_argument.
 *      If the x_i contain duplicates the linear system is ill-conditioned (symbolically it requires a Hermite type).
 *
 * Relation to standard padeApprox: with all x_i = 0 (in theory) it degenerates to the single-point version, but
 * this routine assumes distinct x_i (if you need the x_i = 0 Taylor-coefficient version, use padeApprox).
 *
 * Reference: Stoer & Bulirsch "Introduction to Numerical Analysis" §2.2.4 (rational interp.)
 */
template<concepts::OrderedField T>
[[nodiscard]] RationalInterpResult<T> rationalInterpolation(
    const std::vector<T>& x,
    const std::vector<T>& y,
    std::size_t M, std::size_t N)
{
    using SizeT = typename Matrix<T>::size_type;
    if (x.size() != y.size())
        throw std::invalid_argument("rationalInterpolation: x and y must have same size");
    if (x.size() != M + N + 1)
        throw std::invalid_argument("rationalInterpolation: need exactly M+N+1 data points");

    const std::size_t L = M + N;

    // Build the linear system
    // Variable order: p_0, p_1, ..., p_M, q_1, q_2, ..., q_N
    // Row i: Σ_{k=0..M} p_k x_i^k + Σ_{k=1..N} (-y_i x_i^k) q_k = y_i
    Matrix<T> A(static_cast<SizeT>(L + 1), static_cast<SizeT>(L + 1));
    Vector<T> b(static_cast<SizeT>(L + 1));

    for (std::size_t i = 0; i <= L; ++i) {
        T xp = T(1);
        for (std::size_t k = 0; k <= M; ++k) {
            A(static_cast<SizeT>(i), static_cast<SizeT>(k)) = xp;
            xp = xp * x[i];
        }
        T xq = x[i];  // x_i^1
        for (std::size_t k = 0; k < N; ++k) {
            A(static_cast<SizeT>(i), static_cast<SizeT>(M + 1 + k)) = -y[i] * xq;
            xq = xq * x[i];
        }
        b[static_cast<SizeT>(i)] = y[i];
    }

    Vector<T> sol(static_cast<SizeT>(0));
    try {
        sol = algorithms::solve(A, b, algorithms::SolverType::LU);
    } catch (...) {
        return {{}, {}, false};
    }

    RationalInterpResult<T> r;
    r.numerator.resize(M + 1);
    for (std::size_t k = 0; k <= M; ++k)
        r.numerator[k] = sol[static_cast<SizeT>(k)];
    r.denominator.resize(N + 1);
    r.denominator[0] = T(1);
    for (std::size_t k = 0; k < N; ++k)
        r.denominator[k + 1] = sol[static_cast<SizeT>(M + 1 + k)];

    // Verification: f(x_i) = y_i (within tolerance) at every data point
    // In floating point this is not exact, so verify externally.
    r.valid = true;
    return r;
}

/**
 * @brief Evaluate a rational interpolation result (Horner's method)
 */
template<concepts::OrderedField T>
[[nodiscard]] T evalRational(const RationalInterpResult<T>& r, T x) {
    T num = r.numerator.back();
    for (std::size_t i = r.numerator.size() - 1; i > 0; --i)
        num = num * x + r.numerator[i - 1];
    T den = r.denominator.back();
    for (std::size_t i = r.denominator.size() - 1; i > 0; --i)
        den = den * x + r.denominator[i - 1];
    return num / den;
}

// =====================================================================
// Matrix-valued Padé — matrix exponential e^A (Scaling-and-Squaring + diagonal Padé)
// =====================================================================

/**
 * @brief Compute the matrix exponential e^A.
 *
 * Thin wrapper over the library's production matrix exponential
 * algorithms::expm (scaling-and-squaring with a [13/13] Pade core, plus a
 * Schur-Parlett path for higher precision; Higham 2005). The former stand-alone
 * Pade implementation that lived here was removed in favour of a single,
 * better-tested source of truth.
 *
 * @param A input matrix (square)
 * @param m retained only for source compatibility; the production routine
 *          selects the Pade degree internally, so this argument is ignored.
 * @return e^A
 */
template<typename T>
[[nodiscard]] Matrix<T> matrixExpPade(const Matrix<T>& A, [[maybe_unused]] int m = 6) {
    if (A.rows() != A.cols())
        throw std::invalid_argument("matrixExpPade: matrix must be square");
    return algorithms::expm(A);
}

// =====================================================================
// Least-squares polynomial approximation
// =====================================================================

/**
 * @brief Result of a least-squares polynomial approximation
 */
template<concepts::OrderedField T>
struct PolynomialFitResult {
    std::vector<T> coefficients;   // polynomial coefficients [a0, a1, ..., aN] (a0 + a1*x + ... + aN*x^N)
    T residual;                    // residual sum of squares
};

/**
 * @brief Least-squares polynomial approximation
 *
 * Fit a degree-`degree` polynomial p(x) = a0 + a1*x + ... + a_d*x^d to the data points (x_i, y_i)
 * by least squares.
 *
 * @param x array of independent variables
 * @param y array of dependent variables (same size as x)
 * @param degree polynomial degree (≥ 0)
 * @return PolynomialFitResult{coefficients, residual}
 *
 * @note Improvements:
 *   - mean removal applied for numerical stability (improves the condition number)
 *   - solve the normal equations with algorithms::solve(LU)
 *   - returns the residual sum of squares
 */
template<concepts::OrderedField T>
[[nodiscard]] PolynomialFitResult<T> polynomialFit(
    std::span<const T> x, std::span<const T> y, std::size_t degree) {

    if (x.size() != y.size()) {
        throw std::invalid_argument("polynomialFit: x and y must have the same size");
    }
    if (x.size() < degree + 1) {
        throw std::invalid_argument("polynomialFit: need at least degree+1 data points");
    }

    const std::size_t n = x.size();
    const std::size_t N = degree + 1;  // number of coefficients

    // Subtract the mean for numerical stability
    T mu_x = T(0), mu_y = T(0);
    for (std::size_t i = 0; i < n; ++i) {
        mu_x += x[i];
        mu_y += y[i];
    }
    mu_x /= static_cast<T>(n);
    mu_y /= static_cast<T>(n);

    std::vector<T> xs(n), ys(n);
    for (std::size_t i = 0; i < n; ++i) {
        xs[i] = x[i] - mu_x;
        ys[i] = y[i] - mu_y;
    }

    // Build the normal equations A*c = b
    // A(i,j) = Σ_k xs[k]^(i+j),  b(i) = Σ_k xs[k]^i * ys[k]
    Matrix<T> A(static_cast<typename Matrix<T>::size_type>(N),
                 static_cast<typename Matrix<T>::size_type>(N));
    Vector<T> b(static_cast<typename Vector<T>::size_type>(N));

    // Precompute the powers of xs[k] (up to degree 2*degree)
    std::vector<T> x_pow_sum(2 * degree + 1, T(0));
    for (std::size_t k = 0; k < n; ++k) {
        T xp = T(1);
        for (std::size_t p = 0; p <= 2 * degree; ++p) {
            x_pow_sum[p] += xp;
            xp *= xs[k];
        }
    }

    for (std::size_t i = 0; i < N; ++i) {
        for (std::size_t j = 0; j < N; ++j) {
            A(static_cast<typename Matrix<T>::size_type>(i),
              static_cast<typename Matrix<T>::size_type>(j)) = x_pow_sum[i + j];
        }
        T sum = T(0);
        for (std::size_t k = 0; k < n; ++k) {
            T xp = T(1);
            for (std::size_t p = 0; p < i; ++p) xp *= xs[k];
            sum += xp * ys[k];
        }
        b[static_cast<typename Vector<T>::size_type>(i)] = sum;
    }

    // Solve (coefficients in the shifted space)
    Vector<T> c_shifted = algorithms::solve(A, b, algorithms::SolverType::LU);

    // Convert the shifted-space polynomial p_s(u) = Σ c_k * u^k, u = x - mu_x
    // to the original space p(x) = Σ a_j * x^j
    // p(x) = p_s(x - mu_x) + mu_y
    // (x - mu_x)^k = Σ_{j=0}^{k} C(k,j) * x^j * (-mu_x)^(k-j)

    // Compute the original coefficients by binomial expansion
    std::vector<T> coeffs(N, T(0));
    for (std::size_t k = 0; k < N; ++k) {
        T ck = c_shifted[static_cast<typename Vector<T>::size_type>(k)];
        // Expansion of (x - mu_x)^k
        T binom = T(1);  // C(k, j)
        T neg_mu_pow = T(1);  // (-mu_x)^(k-j) — start at j=0 and update
        // At j=0: C(k,0)*x^0*(-mu_x)^k
        // First compute (-mu_x)^k
        T neg_mu_k = T(1);
        for (std::size_t p = 0; p < k; ++p) neg_mu_k *= (-mu_x);

        for (std::size_t j = 0; j <= k; ++j) {
            // C(k,j) * (-mu_x)^{k-j}
            coeffs[j] += ck * binom * neg_mu_k;
            // Update for the next j: C(k,j+1) = C(k,j) * (k-j)/(j+1)
            //                  (-mu_x)^{k-j-1} = (-mu_x)^{k-j} / (-mu_x)
            if (j < k) {
                binom *= static_cast<T>(k - j) / static_cast<T>(j + 1);
                if (mu_x != T(0)) {
                    neg_mu_k /= (-mu_x);
                } else {
                    neg_mu_k = (k - j - 1 == 0) ? T(1) : T(0);
                }
            }
        }
    }
    coeffs[0] += mu_y;  // add mu_y to the constant term

    // residual sum of squares
    T residual = T(0);
    for (std::size_t i = 0; i < n; ++i) {
        // Evaluate p(x[i]) by Horner's method
        T val = coeffs[degree];
        for (std::size_t j = degree; j > 0; --j) {
            val = val * x[i] + coeffs[j - 1];
        }
        T diff = y[i] - val;
        residual += diff * diff;
    }

    return {coeffs, residual};
}

template<concepts::OrderedField T>
[[nodiscard]] PolynomialFitResult<T> polynomialFit(
    const std::vector<T>& x, const std::vector<T>& y, std::size_t degree) {
    return polynomialFit(std::span<const T>(x), std::span<const T>(y), degree);
}

} // namespace sangi

#endif // SANGI_APPROXIMATION_HPP