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

// Interpolation.hpp
#ifndef INTERPOLATION_HPP
#define INTERPOLATION_HPP

#include <math/concepts/algebraic_concepts.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/linalg/solvers.hpp>
#include <vector>
#include <stdexcept>
#include <array>
#include <span>
#include <optional>

namespace sangi {

/**
 * @brief Evaluate the Lagrange interpolation polynomial
 * @tparam T Field type (e.g. real numbers)
 * @param x Array of x-coordinates of the interpolation points
 * @param y Array of y-coordinates of the interpolation points
 * @param xi Evaluation point
 * @return Interpolated value
 */
template<concepts::Field T>
T lagrange_interpolation(
    std::span<const T> x,
    std::span<const T> y,
    T xi) {

    if (x.size() != y.size() || x.empty()) {
        throw std::invalid_argument("Input arrays must have the same non-zero size");
    }

    T result = T(0);

    for (size_t i = 0; i < x.size(); ++i) {
        T term = y[i];

        for (size_t j = 0; j < x.size(); ++j) {
            if (j != i) {
                term *= (xi - x[j]) / (x[i] - x[j]);
            }
        }

        result += term;
    }

    return result;
}

// Vector version of the Lagrange interpolation polynomial (for backward compatibility)
template<concepts::Field T>
T lagrange_interpolation(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T xi) {
    return lagrange_interpolation(std::span<const T>(x), std::span<const T>(y), xi);
}

/**
 * @brief Compute the coefficients of the Newton interpolation polynomial
 * @tparam T Field type (e.g. real numbers)
 * @param x Array of x-coordinates of the interpolation points
 * @param y Array of y-coordinates of the interpolation points
 * @return Array of divided differences
 */
template<concepts::Field T>
std::vector<T> newton_interpolation_coefficients(
    std::span<const T> x,
    std::span<const T> y) {

    if (x.size() != y.size() || x.empty()) {
        throw std::invalid_argument("Input arrays must have the same non-zero size");
    }

    std::size_t n = x.size();
    std::vector<T> coeffs(y.begin(), y.end());  // Initially the original y values

    // Compute divided differences
    for (std::size_t j = 1; j < n; ++j) {
        for (std::size_t i = n - 1; i >= j; --i) {
            coeffs[i] = (coeffs[i] - coeffs[i - 1]) / (x[i] - x[i - j]);
        }
    }

    return coeffs;
}

// Vector version of the Newton interpolation polynomial coefficient computation (for backward compatibility)
template<concepts::Field T>
std::vector<T> newton_interpolation_coefficients(
    const std::vector<T>& x,
    const std::vector<T>& y) {
    return newton_interpolation_coefficients(std::span<const T>(x), std::span<const T>(y));
}

/**
 * @brief Evaluate the Newton interpolation polynomial
 * @tparam T Field type (e.g. real numbers)
 * @param x Array of x-coordinates of the interpolation points
 * @param coeffs Array of divided-difference coefficients
 * @param xi Evaluation point
 * @return Interpolated value
 */
template<concepts::Field T>
T newton_interpolation_evaluate(
    std::span<const T> x,
    std::span<const T> coeffs,
    T xi) {

    if (x.size() != coeffs.size() || x.empty()) {
        throw std::invalid_argument("Input arrays must have the same non-zero size");
    }

    // Evaluate the polynomial using Horner's method
    T result = coeffs.back();

    for (std::size_t i = coeffs.size() - 1; i > 0; --i) {
        result = result * (xi - x[i - 1]) + coeffs[i - 1];
    }

    return result;
}

// Vector version of the Newton interpolation polynomial evaluation (for backward compatibility)
template<concepts::Field T>
T newton_interpolation_evaluate(
    const std::vector<T>& x,
    const std::vector<T>& coeffs,
    T xi) {
    return newton_interpolation_evaluate(std::span<const T>(x), std::span<const T>(coeffs), xi);
}

// Boundary conditions for the cubic spline
enum class SplineBoundaryCondition {
    Natural,   // Natural spline (second derivative is 0)
    Clamped    // Clamped spline (first-derivative values are specified)
};

/**
 * @brief Compute the coefficients of a cubic spline interpolation
 * @tparam T Field type (e.g. real numbers)
 * @param x Array of x-coordinates of the interpolation points
 * @param y Array of y-coordinates of the interpolation points
 * @param boundary_condition Boundary condition
 * @param bc_values Boundary-condition values (derivative values for a clamped spline)
 * @return Array of spline coefficients
 */
template<concepts::Field T>
std::vector<std::array<T, 4>> cubic_spline_coefficients(
    std::span<const T> x,
    std::span<const T> y,
    SplineBoundaryCondition boundary_condition = SplineBoundaryCondition::Natural,
    const std::pair<T, T>& bc_values = { 0, 0 }) {

    if (x.size() != y.size() || x.size() < 2) {
        throw std::invalid_argument("Input arrays must have the same size >= 2");
    }

    std::size_t n = x.size();
    std::vector<T> h(n - 1);      // Width of each interval
    std::vector<T> alpha(n - 1);  // Temporary right-hand-side values
    std::vector<T> l(n);          // Diagonal entries of the tridiagonal matrix
    std::vector<T> mu(n - 1);     // Super-diagonal entries of the tridiagonal matrix
    std::vector<T> z(n);          // Intermediate solution
    std::vector<T> c(n);          // Spline coefficient c
    std::vector<T> b(n - 1);      // Spline coefficient b
    std::vector<T> d(n - 1);      // Spline coefficient d

    // Compute the step sizes
    for (std::size_t i = 0; i < n - 1; ++i) {
        h[i] = x[i + 1] - x[i];
        if (h[i] <= numeric_traits<T>::zero()) {
            throw std::invalid_argument("X values must be in strictly increasing order");
        }
    }

    // Initialize according to the boundary condition
    if (boundary_condition == SplineBoundaryCondition::Natural) {
        // Natural spline: S''(x0) = S''(xn) = 0
        l[0] = 1;
        mu[0] = 0;
        z[0] = 0;
        l[n - 1] = 1;
        z[n - 1] = 0;
    }
    else {
        // Clamped spline: S'(x0) = f'(x0), S'(xn) = f'(xn)
        l[0] = 2 * h[0];
        mu[0] = 0.5;
        z[0] = 3 * ((y[1] - y[0]) / h[0] - bc_values.first) / h[0];
        l[n - 1] = 2 * h[n - 2];
        z[n - 1] = 3 * (bc_values.second - (y[n - 1] - y[n - 2]) / h[n - 2]) / h[n - 2];
    }

    // Forward elimination of the tridiagonal matrix algorithm
    for (std::size_t i = 1; i < n - 1; ++i) {
        l[i] = 2 * (h[i - 1] + h[i]) - h[i - 1] * mu[i - 1];
        mu[i] = h[i] / l[i];
        alpha[i] = 3 * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]);
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
    }

    // Back substitution
    c[n - 1] = z[n - 1];
    for (std::size_t j = 0; j < n - 1; ++j) {
        std::size_t i = n - 2 - j;
        c[i] = z[i] - mu[i] * c[i + 1];
        b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + 2 * c[i]) / 3;
        d[i] = (c[i + 1] - c[i]) / (3 * h[i]);
    }

    // Store the results in the coefficient array
    std::vector<std::array<T, 4>> coefficients(n - 1);
    for (std::size_t i = 0; i < n - 1; ++i) {
        coefficients[i] = { y[i], b[i], c[i], d[i] };
    }

    return coefficients;
}

// Vector version of the cubic spline interpolation coefficient computation (for backward compatibility)
template<concepts::Field T>
std::vector<std::array<T, 4>> cubic_spline_coefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    int boundary_condition = 0,
    const std::pair<T, T>& bc_values = { 0, 0 }) {
    
    SplineBoundaryCondition bc = (boundary_condition == 0) ? 
        SplineBoundaryCondition::Natural : SplineBoundaryCondition::Clamped;
    
    return cubic_spline_coefficients(
        std::span<const T>(x), 
        std::span<const T>(y), 
        bc, 
        bc_values);
}

/**
 * @brief Evaluate a cubic spline interpolation
 * @tparam T Field type (e.g. real numbers)
 * @param x Array of x-coordinates of the interpolation points
 * @param coeffs Array of spline coefficients
 * @param xi Evaluation point
 * @return Interpolated value
 */
template<concepts::Field T>
T cubic_spline_evaluate(
    std::span<const T> x,
    std::span<const std::array<T, 4>> coeffs,
    T xi) {

    if (x.size() < 2 || coeffs.size() != x.size() - 1) {
        throw std::invalid_argument("Invalid input data for cubic spline evaluation");
    }

    // Search for the interval
    std::size_t i = 0;
    while (i < x.size() - 1 && xi > x[i + 1]) {
        ++i;
    }

    // Out of range: use the value at the endpoint
    if (i >= coeffs.size()) {
        return coeffs.back()[0];  // The a coefficient of the last interval
    }

    // Evaluate the spline
    T dx = xi - x[i];
    const auto& coeff = coeffs[i];
    
    // S(x) = a + b(x-x_i) + c(x-x_i)^2 + d(x-x_i)^3
    return coeff[0] + coeff[1] * dx + coeff[2] * dx * dx + coeff[3] * dx * dx * dx;
}

// Vector version of the cubic spline interpolation evaluation (for backward compatibility)
template<concepts::Field T>
T cubic_spline_evaluate(
    const std::vector<T>& x,
    const std::vector<std::array<T, 4>>& coeffs,
    T xi) {
    return cubic_spline_evaluate(
        std::span<const T>(x), 
        std::span<const std::array<T, 4>>(coeffs), 
        xi);
}

/**
 * @brief Linear interpolation
 * @tparam T Field type (e.g. real numbers)
 * @param x0 x-coordinate 1
 * @param y0 y-coordinate 1
 * @param x1 x-coordinate 2
 * @param y1 y-coordinate 2
 * @param xi Evaluation point
 * @return Interpolated value
 */
template<concepts::Field T>
T linear_interpolate(T x0, T y0, T x1, T y1, T xi) {
    if (approximately_equal<T>(x1, x0)) {
        return (y0 + y1) / 2;  // When x0 == x1, return the average
    }

    T t = (xi - x0) / (x1 - x0);
    return y0 + t * (y1 - y0);
}

/**
 * @brief Bilinear interpolation
 * @tparam T Field type (e.g. real numbers)
 * @param x x-coordinate of the evaluation point
 * @param y y-coordinate of the evaluation point
 * @param x1 x-coordinate of the bottom-left vertex
 * @param y1 y-coordinate of the bottom-left vertex
 * @param x2 x-coordinate of the top-right vertex
 * @param y2 y-coordinate of the top-right vertex
 * @param q11 Value at the bottom-left vertex
 * @param q12 Value at the top-left vertex
 * @param q21 Value at the bottom-right vertex
 * @param q22 Value at the top-right vertex
 * @return Interpolated value
 */
template<concepts::Field T>
T bilinear_interpolate(
    T x, T y,
    T x1, T y1, T x2, T y2,
    T q11, T q12, T q21, T q22) {

    // Check that the coordinates are within range
    if (x < x1 || x > x2 || y < y1 || y > y2) {
        throw std::invalid_argument("Evaluation point outside the interpolation grid");
    }

    T x_ratio = (x - x1) / (x2 - x1);
    T y_ratio = (y - y1) / (y2 - y1);

    // Linear interpolation along the bottom and top edges
    T r1 = linear_interpolate(T(0), q11, T(1), q21, x_ratio);
    T r2 = linear_interpolate(T(0), q12, T(1), q22, x_ratio);

    // Linear interpolation of the results in the vertical direction
    return linear_interpolate(T(0), r1, T(1), r2, y_ratio);
}

// =====================================================================
// B-spline interpolation
// =====================================================================

/**
 * @brief Recursive evaluation of the B-spline basis function (Cox-de Boor)
 * @param degree Degree of the B-spline
 * @param i Index of the basis function
 * @param x Evaluation point
 * @param knots Knot vector
 * @return Value of B_{i,degree}(x)
 */
template<concepts::OrderedField T>
[[nodiscard]] T bsplineBasis(
    std::size_t degree,
    std::size_t i,
    T x,
    std::span<const T> knots) {

    // Degree-0 basis
    if (degree == 0) {
        // The last knot interval includes the right endpoint
        if (i + 1 < knots.size() && knots[i + 1] == knots.back()) {
            return (knots[i] <= x && x <= knots[i + 1]) ? T(1) : T(0);
        }
        return (knots[i] <= x && x < knots[i + 1]) ? T(1) : T(0);
    }

    T left = T(0), right = T(0);

    T denom_left = knots[i + degree] - knots[i];
    if (denom_left > T(0)) {
        left = (x - knots[i]) / denom_left * bsplineBasis(degree - 1, i, x, knots);
    }

    T denom_right = knots[i + degree + 1] - knots[i + 1];
    if (denom_right > T(0)) {
        right = (knots[i + degree + 1] - x) / denom_right * bsplineBasis(degree - 1, i + 1, x, knots);
    }

    return left + right;
}

/**
 * @brief Compute the knot vector and control-point coefficients for B-spline interpolation
 * @param x x-coordinates of the data points (ascending, n of them)
 * @param y y-coordinates of the data points (same size as x)
 * @param degree Degree of the B-spline (default 3: cubic B-spline)
 * @return The pair {knots, coefficients}
 * @throw std::invalid_argument Insufficient data or size mismatch
 *
 * @note Uses a clamped knot vector:
 *   - Place knots with multiplicity (degree+1) at both ends
 *   - Interior knots use averaged-knot selection
 *
 * @note Improvements:
 *   - Set the boundary rows explicitly (to prevent numerical error)
 *   - Use algorithms::solve(LU)
 */
template<concepts::OrderedField T>
[[nodiscard]] std::pair<std::vector<T>, std::vector<T>> bsplineCoefficients(
    std::span<const T> x,
    std::span<const T> y,
    std::size_t degree = 3) {

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

    // Build the clamped knot vector
    // Total number of knots = n + degree + 1
    std::size_t num_knots = n + degree + 1;
    std::vector<T> knots(num_knots);

    // Beginning: multiplicity (degree+1)
    for (std::size_t i = 0; i <= degree; ++i) {
        knots[i] = x.front();
    }
    // End: multiplicity (degree+1)
    for (std::size_t i = 0; i <= degree; ++i) {
        knots[num_knots - 1 - i] = x.back();
    }
    // Interior knots: averaged knots
    for (std::size_t j = 1; j < n - degree; ++j) {
        T sum = T(0);
        for (std::size_t i = j; i < j + degree; ++i) {
            sum += x[i];
        }
        knots[j + degree] = sum / static_cast<T>(degree);
    }

    // Build the collocation matrix
    // The n×n matrix of B_{j,degree}(x_i)
    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));

    std::span<const T> knots_span(knots);

    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)) = bsplineBasis(degree, j, x[i], knots_span);
        }
        b[static_cast<typename Vector<T>::size_type>(i)] = y[i];
    }

    // Set the boundary rows explicitly (to prevent numerical error)
    // At x[0], B_{0,degree}(x[0]) = 1, all others are 0
    for (std::size_t j = 0; j < n; ++j) {
        A(0, static_cast<typename Matrix<T>::size_type>(j)) = T(0);
    }
    A(0, 0) = T(1);
    b[0] = y[0];

    // At x[n-1], B_{n-1,degree}(x[n-1]) = 1, all others are 0
    auto last = static_cast<typename Matrix<T>::size_type>(n - 1);
    for (std::size_t j = 0; j < n; ++j) {
        A(last, static_cast<typename Matrix<T>::size_type>(j)) = T(0);
    }
    A(last, last) = T(1);
    b[static_cast<typename Vector<T>::size_type>(n - 1)] = y[n - 1];

    // Solve the linear system
    Vector<T> coeffs = algorithms::solve(A, b, algorithms::SolverType::LU);

    std::vector<T> result(n);
    for (std::size_t i = 0; i < n; ++i) {
        result[i] = coeffs[static_cast<typename Vector<T>::size_type>(i)];
    }

    return {knots, result};
}

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

/**
 * @brief Evaluate a B-spline interpolation
 * @param knots Knot vector
 * @param coefficients Control-point coefficients
 * @param degree Degree of the B-spline
 * @param x Evaluation point
 * @return Interpolated value
 */
template<concepts::OrderedField T>
[[nodiscard]] T bsplineEvaluate(
    std::span<const T> knots,
    std::span<const T> coefficients,
    std::size_t degree,
    T x) {

    T result = T(0);
    for (std::size_t i = 0; i < coefficients.size(); ++i) {
        result += coefficients[i] * bsplineBasis(degree, i, x, knots);
    }
    return result;
}

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

// =====================================================================
// Barycentric Lagrange interpolation
// =====================================================================

/**
 * @brief Compute the weights for barycentric Lagrange interpolation
 * @param x Array of x-coordinates of the interpolation points (elements must be distinct)
 * @return Barycentric weights w[j] = 1 / Π_{i≠j}(x[j] - x[i])
 *
 * @note Ordinary Lagrange interpolation is O(n^2) every time, but the barycentric form is
 *       O(n^2) for the weight computation + O(n) per evaluation. Advantageous when evaluating
 *       multiple times with the same interpolation points.
 */
template<concepts::Field T>
[[nodiscard]] std::vector<T> barycentricLagrangeWeights(std::span<const T> x) {
    if (x.empty()) {
        throw std::invalid_argument("barycentricLagrangeWeights: empty input");
    }

    std::size_t n = x.size();
    std::vector<T> w(n);

    for (std::size_t j = 0; j < n; ++j) {
        T wj = T(1);
        for (std::size_t i = 0; i < n; ++i) {
            if (i != j) {
                wj *= (x[j] - x[i]);
            }
        }
        w[j] = T(1) / wj;
    }

    return w;
}

template<concepts::Field T>
[[nodiscard]] std::vector<T> barycentricLagrangeWeights(const std::vector<T>& x) {
    return barycentricLagrangeWeights(std::span<const T>(x));
}

/**
 * @brief Evaluate barycentric Lagrange interpolation
 * @param x Array of x-coordinates of the interpolation points
 * @param y Array of y-coordinates of the interpolation points
 * @param weights Weights computed by barycentricLagrangeWeights()
 * @param xi Evaluation point
 * @return Interpolated value
 */
template<concepts::Field T>
[[nodiscard]] T barycentricLagrangeEvaluate(
    std::span<const T> x,
    std::span<const T> y,
    std::span<const T> weights,
    T xi) {

    if (x.size() != y.size() || x.size() != weights.size() || x.empty()) {
        throw std::invalid_argument("barycentricLagrangeEvaluate: size mismatch or empty input");
    }

    std::size_t n = x.size();
    T num = T(0);
    T den = T(0);

    for (std::size_t j = 0; j < n; ++j) {
        T d = xi - x[j];
        if (d == T(0)) return y[j];  // Coincides with an interpolation point

        T c = weights[j] / d;
        num += c * y[j];
        den += c;
    }

    return num / den;
}

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

// =====================================================================
// Step interpolation (Zero-order hold)
// =====================================================================

/**
 * @brief Step interpolation (zero-order hold)
 * @param x Array of x-coordinates of the interpolation points (ascending)
 * @param y Array of y-coordinates of the interpolation points
 * @param xi Evaluation point
 * @return y[i] at the left end of the interval [x[i], x[i+1]) containing xi
 *
 * @note Out of range: returns the value at the nearest endpoint.
 */
template<concepts::OrderedField T>
[[nodiscard]] T stepInterpolation(
    std::span<const T> x,
    std::span<const T> y,
    T xi) {

    if (x.size() != y.size() || x.empty()) {
        throw std::invalid_argument("stepInterpolation: size mismatch or empty input");
    }

    std::size_t n = x.size();

    if (n == 1) return y[0];

    // Out of range
    if (xi <= x[0]) return y[0];
    if (xi >= x[n - 1]) return y[n - 1];

    // Binary search for the interval: x[lo] <= xi < x[lo+1]
    std::size_t lo = 0, hi = n - 1;
    while (lo + 1 < hi) {
        std::size_t mid = (lo + hi) / 2;
        if (x[mid] <= xi) lo = mid;
        else hi = mid;
    }

    return y[lo];
}

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

// =====================================================================
// Smooth piecewise interpolation
// =====================================================================

/**
 * @brief Compute the coefficients of smooth interpolation
 *
 * Interpolates each interval with a piecewise polynomial. A Hermite-type
 * interpolation that mixes the slopes of adjacent intervals to achieve a smooth join.
 * - First interval: quadratic (constrains the slope at the end point)
 * - Middle intervals: cubic (constrains the slopes at both ends)
 * - Last interval: quadratic (constrains the slope at the start point)
 *
 * Slope mixing formula: Amix = (|A0|·A1 + |A1|·A0) / (|A0|+|A1|)
 * - 3 points on a straight line -> reproduces the straight line
 * - One side has slope 0 (flat interval) -> the mix is also 0
 * - Slopes have opposite signs (peak) -> the mix is 0
 *
 * @param x Array of x-coordinates of the interpolation points (ascending, 2 or more points)
 * @param y Array of y-coordinates of the interpolation points
 * @return Coefficients of each interval {c0, c1, c2, c3}: f(t) = c0 + c1·t + c2·t² + c3·t³, t = x - x[i]
 *
 * @note cubic_spline_evaluate() can be used for evaluation (the coefficient form is identical).
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 4>> smoothInterpolationCoefficients(
    std::span<const T> x,
    std::span<const T> y) {

    if (x.size() != y.size() || x.size() < 2) {
        throw std::invalid_argument(
            "smoothInterpolationCoefficients: need at least 2 points with matching sizes");
    }

    std::size_t n = x.size();
    std::size_t intervals = n - 1;

    // Slope of each interval
    std::vector<T> slope(intervals);
    for (std::size_t i = 0; i < intervals; ++i) {
        T dx = x[i + 1] - x[i];
        if (dx <= T(0)) {
            throw std::invalid_argument(
                "smoothInterpolationCoefficients: x must be strictly increasing");
        }
        slope[i] = (y[i + 1] - y[i]) / dx;
    }

    // Slopes at the start/end points of each interval (initial value is the interval slope)
    std::vector<T> A0(intervals), A1(intervals);
    for (std::size_t i = 0; i < intervals; ++i) {
        A0[i] = A1[i] = slope[i];
    }

    // Mix the slopes at the interior nodes
    for (std::size_t i = 1; i < intervals; ++i) {
        T a0 = slope[i];      // Slope of the current interval
        T a1 = slope[i - 1];  // Slope of the previous interval
        T abs_a0 = (a0 < T(0)) ? -a0 : a0;
        T abs_a1 = (a1 < T(0)) ? -a1 : a1;
        T sum_abs = abs_a0 + abs_a1;
        T mix = T(0);
        if (sum_abs > T(0)) {
            mix = (abs_a0 * a1 + abs_a1 * a0) / sum_abs;
        }
        A1[i - 1] = mix;
        A0[i] = mix;
    }

    // Build the polynomial coefficients for each interval
    std::vector<std::array<T, 4>> coeffs(intervals);
    for (std::size_t i = 0; i < intervals; ++i) {
        T y0 = y[i];
        T y1 = y[i + 1];
        T x1 = x[i + 1] - x[i];
        T a0 = A0[i];
        T a1 = A1[i];
        T x12 = x1 * x1;

        if (i == 0 && intervals > 1) {
            // First interval: quadratic (constrains the end-point slope a1)
            T ca = (a1 * x1 + y0 - y1) / x12;
            T cb = -(a1 * x12 + T(2) * x1 * (y0 - y1)) / x12;
            coeffs[i] = { y0, cb, ca, T(0) };
        } else if (i == intervals - 1 && intervals > 1) {
            // Last interval: quadratic (constrains the start-point slope a0)
            T ca = -(a0 * x1 + y0 - y1) / x12;
            coeffs[i] = { y0, a0, ca, T(0) };
        } else {
            // Middle interval (or single interval): cubic (constrains the slopes at both ends)
            T x13 = x12 * x1;
            T ca = ((a0 + a1) * x1 + T(2) * (y0 - y1)) / x13;
            T cb = -((T(2) * a0 + a1) * x12 + T(3) * x1 * (y0 - y1)) / x13;
            coeffs[i] = { y0, a0, cb, ca };
        }
    }

    return coeffs;
}

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

// =====================================================================
// PCHIP (Piecewise Cubic Hermite Interpolating Polynomial)
// =====================================================================

/**
 * @brief Compute PCHIP interpolation coefficients
 *
 * Monotonicity-preserving interpolation via the Fritsch-Carlson algorithm.
 * Preserves the shape of rising/falling edges and suppresses spurious oscillation.
 * The return value can be evaluated with cubic_spline_evaluate().
 *
 * @param x x-coordinates of the nodes (ascending, n≥2)
 * @param y y-coordinates of the nodes
 * @return Polynomial coefficients {a, b, c, d} for each interval
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 4>> pchipCoefficients(
    std::span<const T> x, std::span<const T> y)
{
    std::size_t n = x.size();
    if (n < 2 || x.size() != y.size()) return {};

    std::size_t m = n - 1;  // Number of intervals

    // Interval widths and slopes
    std::vector<T> h(m), delta(m);
    for (std::size_t i = 0; i < m; ++i) {
        h[i] = x[i + 1] - x[i];
        delta[i] = (y[i + 1] - y[i]) / h[i];
    }

    // Slope d[k] at each node
    std::vector<T> d(n, T{0});

    if (n == 2) {
        d[0] = delta[0];
        d[1] = delta[0];
    } else {
        // Interior points: Fritsch-Carlson harmonic mean
        for (std::size_t i = 1; i < m; ++i) {
            if (delta[i - 1] * delta[i] > T{0}) {
                // Same sign: weighted harmonic mean
                T w1 = T{2} * h[i] + h[i - 1];
                T w2 = h[i] + T{2} * h[i - 1];
                d[i] = (w1 + w2) / (w1 / delta[i - 1] + w2 / delta[i]);
            }
            // else: signs differ or are 0 -> d[i] = 0 (default value)
        }

        // Endpoints: Bessel's one-sided 3-point formula
        d[0] = ((T{2} * h[0] + h[1]) * delta[0] - h[0] * delta[1]) / (h[0] + h[1]);
        // Monotonicity constraint: 0 if d[0] and delta[0] have opposite signs
        if (d[0] * delta[0] < T{0}) {
            d[0] = T{0};
        } else if (delta[0] * delta[1] < T{0} && std::abs(d[0]) > T{3} * std::abs(delta[0])) {
            d[0] = T{3} * delta[0];
        }

        d[m] = ((T{2} * h[m - 1] + h[m - 2]) * delta[m - 1] - h[m - 1] * delta[m - 2]) / (h[m - 1] + h[m - 2]);
        if (d[m] * delta[m - 1] < T{0}) {
            d[m] = T{0};
        } else if (delta[m - 1] * delta[m - 2] < T{0} && std::abs(d[m]) > T{3} * std::abs(delta[m - 1])) {
            d[m] = T{3} * delta[m - 1];
        }

        // Fritsch-Carlson monotonicity condition: α² + β² ≤ 9
        for (std::size_t i = 0; i < m; ++i) {
            if (delta[i] == T{0}) {
                d[i] = T{0};
                d[i + 1] = T{0};
            } else {
                T alpha = d[i] / delta[i];
                T beta = d[i + 1] / delta[i];
                T r2 = alpha * alpha + beta * beta;
                if (r2 > T{9}) {
                    T tau = T{3} / std::sqrt(r2);
                    d[i] = tau * alpha * delta[i];
                    d[i + 1] = tau * beta * delta[i];
                }
            }
        }
    }

    // Convert Hermite -> polynomial coefficients array<T,4>
    // S(x) = a + b(x-xi) + c(x-xi)² + d(x-xi)³
    // Hermite: p(t) = y_i(1-t)²(1+2t) + y_{i+1}t²(3-2t) + d_i·h·t(1-t)² + d_{i+1}·h·t²(t-1)
    // where t = (x-xi)/h
    std::vector<std::array<T, 4>> coeffs(m);
    for (std::size_t i = 0; i < m; ++i) {
        T hi = h[i];
        coeffs[i][0] = y[i];                                     // a = y_i
        coeffs[i][1] = d[i];                                     // b = d_i
        coeffs[i][2] = (T{3} * delta[i] - T{2} * d[i] - d[i + 1]) / hi;    // c
        coeffs[i][3] = (d[i] + d[i + 1] - T{2} * delta[i]) / (hi * hi);    // d
    }

    return coeffs;
}

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

// =====================================================================
// Floater-Hormann barycentric rational interpolation
// =====================================================================

/**
 * @brief Compute the weights for Floater-Hormann barycentric rational interpolation
 *
 * Blended polynomial interpolation of degree d. Avoids the Runge phenomenon while
 * remaining stable even at equally spaced points.
 * d=0: piecewise constant, d=n-1: polynomial interpolation
 *
 * @param x x-coordinates of the nodes (n points)
 * @param d Blend degree (0 ≤ d ≤ n-1)
 * @return Weight w[k] for each node
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<T> floaterHormannWeights(
    std::span<const T> x, std::size_t d = 3)
{
    std::size_t n = x.size();
    if (n == 0) return {};
    if (d >= n) d = n - 1;

    std::vector<T> w(n, T{0});

    for (std::size_t k = 0; k < n; ++k) {
        T sum = T{0};
        std::size_t i_min = (k >= d) ? k - d : 0;
        std::size_t i_max = std::min(k, n - 1 - d);

        for (std::size_t i = i_min; i <= i_max; ++i) {
            T prod = T{1};
            for (std::size_t j = i; j <= i + d; ++j) {
                if (j != k) {
                    prod /= (x[k] - x[j]);
                }
            }
            // Sign of (-1)^i (Floater-Hormann definition: applied to the inner sum)
            T sign_i = (i % 2 == 0) ? T{1} : T{-1};
            sum += sign_i * prod;
        }

        w[k] = sum;
    }

    return w;
}

/**
 * @brief Evaluate Floater-Hormann barycentric rational interpolation
 *
 * Barycentric formula: r(xi) = Σ w[k]·y[k]/(xi-x[k]) / Σ w[k]/(xi-x[k])
 * Returns exact values at the nodes.
 */
template<concepts::OrderedField T>
[[nodiscard]] T floaterHormannEvaluate(
    std::span<const T> x, std::span<const T> y,
    std::span<const T> weights, T xi)
{
    std::size_t n = x.size();

    T num = T{0};
    T den = T{0};

    for (std::size_t k = 0; k < n; ++k) {
        T diff = xi - x[k];
        if (std::abs(diff) < std::numeric_limits<T>::epsilon() * T{10}) {
            return y[k];  // On a node
        }
        T term = weights[k] / diff;
        num += term * y[k];
        den += term;
    }

    return num / den;
}

// vector convenience overload
template<concepts::OrderedField T>
[[nodiscard]] std::vector<T> floaterHormannWeights(
    const std::vector<T>& x, std::size_t d = 3) {
    return floaterHormannWeights(std::span<const T>(x), d);
}

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

// =====================================================================
// Catmull-Rom spline
// =====================================================================

/**
 * @brief Compute Catmull-Rom spline coefficients
 *
 * A C¹-continuous interpolating spline. The tangent at each node is computed from
 * its two neighbors. A standard curve-interpolation method in computer graphics.
 * The return value can be evaluated with cubic_spline_evaluate().
 *
 * @param x x-coordinates of the nodes (ascending, n≥2)
 * @param y y-coordinates of the nodes
 * @return Polynomial coefficients {a, b, c, d} for each interval
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 4>> catmullRomCoefficients(
    std::span<const T> x, std::span<const T> y)
{
    std::size_t n = x.size();
    if (n < 2 || x.size() != y.size()) return {};

    std::size_t m = n - 1;  // Number of intervals

    // Tangent at each node
    std::vector<T> t(n);
    if (n == 2) {
        t[0] = (y[1] - y[0]) / (x[1] - x[0]);
        t[1] = t[0];
    } else {
        // Endpoints: one-sided difference
        t[0] = (y[1] - y[0]) / (x[1] - x[0]);
        t[n - 1] = (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]);

        // Interior points: difference of the two neighbors
        for (std::size_t i = 1; i < n - 1; ++i) {
            t[i] = (y[i + 1] - y[i - 1]) / (x[i + 1] - x[i - 1]);
        }
    }

    // Hermite basis -> polynomial coefficients
    std::vector<std::array<T, 4>> coeffs(m);
    for (std::size_t i = 0; i < m; ++i) {
        T hi = x[i + 1] - x[i];
        T delta = (y[i + 1] - y[i]) / hi;

        coeffs[i][0] = y[i];                                           // a
        coeffs[i][1] = t[i];                                           // b
        coeffs[i][2] = (T{3} * delta - T{2} * t[i] - t[i + 1]) / hi;  // c
        coeffs[i][3] = (t[i] + t[i + 1] - T{2} * delta) / (hi * hi);  // d
    }

    return coeffs;
}

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

// =====================================================================
// Modified Akima interpolation
// =====================================================================

/**
 * @brief Compute Modified Akima interpolation coefficients
 *
 * Applies a Boost.Math-conformant modification to Akima's local slope estimation.
 * Adding ε to the absolute value of the δ difference avoids division by zero and
 * improves stability. The return value can be evaluated with cubic_spline_evaluate().
 *
 * @param x x-coordinates of the nodes (ascending, n≥2)
 * @param y y-coordinates of the nodes
 * @return Polynomial coefficients {a, b, c, d} for each interval
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 4>> modifiedAkimaCoefficients(
    std::span<const T> x, std::span<const T> y)
{
    std::size_t n = x.size();
    if (n < 2 || x.size() != y.size()) return {};

    if (n == 2) {
        T slope = (y[1] - y[0]) / (x[1] - x[0]);
        return { {y[0], slope, T{0}, T{0}} };
    }

    std::size_t m = n - 1;  // Number of intervals

    // δ[i] = (y[i+1]-y[i]) / (x[i+1]-x[i])
    std::vector<T> delta(m);
    for (std::size_t i = 0; i < m; ++i) {
        delta[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
    }

    // Generate virtual δ values for the endpoints by linear extrapolation
    // δ[-2], δ[-1], ..., δ[m], δ[m+1]
    // m+4 in total: idx 0,1 are left-end extrapolation, idx 2..m+1 are the body, idx m+2,m+3 are right-end extrapolation
    std::vector<T> d(m + 4);
    for (std::size_t i = 0; i < m; ++i) {
        d[i + 2] = delta[i];
    }
    // Left-end extrapolation
    d[1] = T{2} * delta[0] - (m > 1 ? delta[1] : delta[0]);
    d[0] = T{2} * d[1] - delta[0];
    // Right-end extrapolation
    d[m + 2] = T{2} * delta[m - 1] - (m > 1 ? delta[m - 2] : delta[m - 1]);
    d[m + 3] = T{2} * d[m + 2] - delta[m - 1];

    // Modified Akima weight computation
    // w_i = |δ_{i+1} - δ_i| + ε  (stabilized by ε)
    T max_abs_delta = T{0};
    for (std::size_t i = 0; i < m + 3; ++i) {
        T ad = std::abs(d[i + 1] - d[i]);
        if (ad > max_abs_delta) max_abs_delta = ad;
    }
    T eps = std::numeric_limits<T>::epsilon() * max_abs_delta;
    if (eps == T{0}) eps = std::numeric_limits<T>::epsilon();

    // Slope t[i] at each node
    std::vector<T> t(n);
    for (std::size_t i = 0; i < n; ++i) {
        // Indices into d: for node i, d[i], d[i+1], d[i+2], d[i+3]
        T w1 = std::abs(d[i + 3] - d[i + 2]) + eps;
        T w2 = std::abs(d[i + 1] - d[i]) + eps;
        t[i] = (w1 * d[i + 1] + w2 * d[i + 2]) / (w1 + w2);
    }

    // Hermite -> polynomial coefficients
    std::vector<std::array<T, 4>> coeffs(m);
    for (std::size_t i = 0; i < m; ++i) {
        T hi = x[i + 1] - x[i];
        T di = delta[i];

        coeffs[i][0] = y[i];
        coeffs[i][1] = t[i];
        coeffs[i][2] = (T{3} * di - T{2} * t[i] - t[i + 1]) / hi;
        coeffs[i][3] = (t[i] + t[i + 1] - T{2} * di) / (hi * hi);
    }

    return coeffs;
}

// vector convenience overload
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 4>> modifiedAkimaCoefficients(
    const std::vector<T>& x, const std::vector<T>& y) {
    return modifiedAkimaCoefficients(std::span<const T>(x), std::span<const T>(y));
}

// =====================================================================
// Quintic Hermite interpolation
// =====================================================================

/**
 * @brief Compute Quintic Hermite interpolation coefficients
 *
 * Quintic Hermite interpolation specifying f, f', f'' at each node. C²-continuous.
 * 6 coefficients {a,b,c,d,e,f}: S(t) = a + b·t + c·t² + d·t³ + e·t⁴ + f·t⁵ (t = x-x_i)
 *
 * @param x   x-coordinates of the nodes (ascending, n≥2)
 * @param y   y-values of the nodes
 * @param dy  First-derivative value at each node
 * @param d2y Second-derivative value at each node
 * @return Polynomial coefficients array<T,6> for each interval
 */
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 6>> quinticHermiteCoefficients(
    std::span<const T> x, std::span<const T> y,
    std::span<const T> dy, std::span<const T> d2y)
{
    std::size_t n = x.size();
    if (n < 2 || y.size() != n || dy.size() != n || d2y.size() != n) return {};

    std::size_t m = n - 1;
    std::vector<std::array<T, 6>> coeffs(m);

    for (std::size_t i = 0; i < m; ++i) {
        T h = x[i + 1] - x[i];
        T h2 = h * h;
        T h3 = h2 * h;
        T h4 = h3 * h;
        T h5 = h4 * h;

        T y0 = y[i], y1 = y[i + 1];
        T dy0 = dy[i], dy1 = dy[i + 1];
        T d2y0 = d2y[i], d2y1 = d2y[i + 1];

        // Derive the coefficients from the quintic Hermite basis
        // S(0) = y0, S(h) = y1
        // S'(0) = dy0, S'(h) = dy1
        // S''(0) = d2y0, S''(h) = d2y1
        coeffs[i][0] = y0;
        coeffs[i][1] = dy0;
        coeffs[i][2] = d2y0 / T{2};

        // a3, a4, a5 come from the system of equations
        // S(h) = y0 + dy0·h + d2y0/2·h² + a3·h³ + a4·h⁴ + a5·h⁵ = y1
        // S'(h) = dy0 + d2y0·h + 3a3·h² + 4a4·h³ + 5a5·h⁴ = dy1
        // S''(h) = d2y0 + 6a3·h + 12a4·h² + 20a5·h³ = d2y1

        T dy_diff = y1 - y0 - dy0 * h - d2y0 / T{2} * h2;
        T ddy_diff = dy1 - dy0 - d2y0 * h;
        T dddy_diff = d2y1 - d2y0;

        // System of equations (Cramer's rule):
        // a3·h³ + a4·h⁴ + a5·h⁵ = dy_diff
        // 3a3·h² + 4a4·h³ + 5a5·h⁴ = ddy_diff
        // 6a3·h + 12a4·h² + 20a5·h³ = dddy_diff
        // Normalize by dividing through by h:
        // a3 + a4·h + a5·h² = dy_diff/h³
        // 3a3 + 4a4·h + 5a5·h² = ddy_diff/h²
        // 6a3 + 12a4·h + 20a5·h² = dddy_diff/h

        T r1 = dy_diff / h3;
        T r2 = ddy_diff / h2;
        T r3 = dddy_diff / h;

        // Solution: Gaussian elimination
        // Row2 - 3*Row1: a4·h + 2a5·h² = r2 - 3r1
        // Row3 - 6*Row1: 6a4·h + 14a5·h² = r3 - 6r1
        T s1 = r2 - T{3} * r1;
        T s2 = r3 - T{6} * r1;
        // Row3' - 6*Row2': 2a5·h² = s2 - 6s1
        T a5 = (s2 - T{6} * s1) / (T{2} * h2);
        T a4 = (s1 - T{2} * a5 * h2) / h;
        T a3 = r1 - a4 * h - a5 * h2;

        coeffs[i][3] = a3;
        coeffs[i][4] = a4;
        coeffs[i][5] = a5;
    }

    return coeffs;
}

// vector convenience overload
template<concepts::OrderedField T>
[[nodiscard]] std::vector<std::array<T, 6>> quinticHermiteCoefficients(
    const std::vector<T>& x, const std::vector<T>& y,
    const std::vector<T>& dy, const std::vector<T>& d2y) {
    return quinticHermiteCoefficients(
        std::span<const T>(x), std::span<const T>(y),
        std::span<const T>(dy), std::span<const T>(d2y));
}

/**
 * @brief Evaluate Quintic Hermite interpolation
 *
 * @param x      x-coordinates of the nodes
 * @param coeffs Return value of quinticHermiteCoefficients()
 * @param xi     Evaluation point
 * @return Interpolated value
 */
template<concepts::OrderedField T>
[[nodiscard]] T quinticHermiteEvaluate(
    std::span<const T> x,
    std::span<const std::array<T, 6>> coeffs,
    T xi)
{
    std::size_t n = x.size();
    if (n < 2 || coeffs.size() != n - 1) return T{0};

    // Binary search for the interval
    std::size_t idx = 0;
    if (xi <= x[0]) {
        idx = 0;
    } else if (xi >= x[n - 1]) {
        idx = n - 2;
    } else {
        auto it = std::upper_bound(x.begin(), x.end(), xi);
        idx = static_cast<std::size_t>(std::distance(x.begin(), it)) - 1;
        if (idx >= n - 1) idx = n - 2;
    }

    T t = xi - x[idx];
    const auto& c = coeffs[idx];
    return c[0] + t * (c[1] + t * (c[2] + t * (c[3] + t * (c[4] + t * c[5]))));
}

template<concepts::OrderedField T>
[[nodiscard]] T quinticHermiteEvaluate(
    const std::vector<T>& x,
    const std::vector<std::array<T, 6>>& coeffs,
    T xi) {
    return quinticHermiteEvaluate(std::span<const T>(x),
        std::span<const std::array<T, 6>>(coeffs), xi);
}

// =====================================================================
// Cardinal Trigonometric Interpolation
// =====================================================================

/**
 * @brief Cardinal trigonometric interpolation
 *
 * Interpolates equally spaced periodic data y[0..N-1] with a trigonometric polynomial.
 * An interpolation method that underpins periodic functions and spectral methods.
 * Direct-computation approach (no FFT, suited to small-to-medium data).
 *
 * @param y      Equally spaced sample values y[0..N-1]
 * @param period Period T (x[k] = k·period/N)
 * @param xi     Evaluation point
 * @return Interpolated value
 */
template<concepts::OrderedField T>
[[nodiscard]] T cardinalTrigonometricInterpolate(
    std::span<const T> y, T period, T xi)
{
    std::size_t N = y.size();
    if (N == 0) return T{0};
    if (N == 1) return y[0];

    const T pi = std::acos(T{-1});
    T sum = T{0};

    // Dirichlet kernel approach:
    // S(x) = (1/N) Σ_{k=0}^{N-1} y[k] · D_N(x - x_k)
    // D_N(t) = sin(N·π·t/P) / (N·sin(π·t/P))  (N odd)
    //        = sin(N·π·t/P) / (N·tan(π·t/P))  (N even, no Fejér-kernel-style correction needed)
    // where D_N=1 when t=0

    for (std::size_t k = 0; k < N; ++k) {
        T xk = static_cast<T>(k) * period / static_cast<T>(N);
        T t = xi - xk;

        // Period normalization
        T arg = pi * t / period;
        T sin_arg = std::sin(arg);

        T dk;
        if (std::abs(sin_arg) < std::numeric_limits<T>::epsilon() * T{100}) {
            // t ≈ 0 (mod period): kernel value = 1
            // Decide by the sign of cos(Nπt/P) (may be -1 when t ≈ P/2)
            T cos_narg = std::cos(static_cast<T>(N) * arg);
            dk = (cos_narg > T{0}) ? T{1} : T{-1};
        } else if (N % 2 == 0) {
            // Even N: φ(t) = sin(Nπt/P) / (N·tan(πt/P))
            T sin_narg = std::sin(static_cast<T>(N) * arg);
            T tan_arg = std::tan(arg);
            dk = sin_narg / (static_cast<T>(N) * tan_arg);
        } else {
            // Odd N: φ(t) = sin(Nπt/P) / (N·sin(πt/P))
            T sin_narg = std::sin(static_cast<T>(N) * arg);
            dk = sin_narg / (static_cast<T>(N) * sin_arg);
        }
        sum += y[k] * dk;
    }

    return sum;
}

// vector convenience overload
template<concepts::OrderedField T>
[[nodiscard]] T cardinalTrigonometricInterpolate(
    const std::vector<T>& y, T period, T xi) {
    return cardinalTrigonometricInterpolate(std::span<const T>(y), period, xi);
}

// MKL implementation (if available)
// =====================================================================
// FITPACK: spline curve/surface fitting (Dierckx method)
// =====================================================================

/**
 * @brief Result of smoothing spline fitting
 */
template<concepts::OrderedField T>
struct SplineFitResult {
    std::vector<T> knots;        ///< Knot vector
    std::vector<T> coefficients; ///< B-spline coefficients
    std::size_t degree;          ///< Spline degree
    T smoothing;                 ///< The actual smoothing parameter
    T residual;                  ///< Sum of squared residuals
};

/**
 * @brief Smoothing spline fitting (equivalent to FITPACK splrep)
 *
 * Performs B-spline fitting with a smoothing condition on the data points (x_i, y_i).
 *
 * Minimize: Σ w_i (y_i - S(x_i))² + s · ∫ (S''(x))² dx
 *   s = 0: interpolation (passes exactly through the data points)
 *   s > 0: smoothing (reduces the effect of noise)
 *
 * The number of knots is selected automatically according to the smoothing parameter.
 *
 * @param x       x-coordinates of the data points (ascending, n of them)
 * @param y       y-coordinates of the data points (n of them)
 * @param w       Weights (all 1.0 if empty)
 * @param s       Smoothing parameter (0 = interpolation, < 0 = automatic)
 * @param degree  Spline degree (default 3 = cubic)
 * @return SplineFitResult<T>
 * @throws std::invalid_argument Insufficient data
 */
template<concepts::OrderedField T>
[[nodiscard]]
SplineFitResult<T> splineFit(
    std::span<const T> x,
    std::span<const T> y,
    std::span<const T> w = {},
    T s = T(-1),
    std::size_t degree = 3)
{
    std::size_t n = x.size();
    if (n < degree + 1)
        throw std::invalid_argument("splineFit: need at least (degree+1) data points");
    if (y.size() != n)
        throw std::invalid_argument("splineFit: x and y must have same size");
    if (!w.empty() && w.size() != n)
        throw std::invalid_argument("splineFit: w must be empty or same size as x");

    // Set up the weights
    std::vector<T> weights(n, T(1));
    if (!w.empty()) {
        for (std::size_t i = 0; i < n; ++i) weights[i] = w[i];
    }

    // Automatic smoothing parameter: s = n (Dierckx's recommendation)
    if (s < T(0)) s = static_cast<T>(n);

    // When s == 0: interpolation (use the existing bsplineCoefficients)
    if (s == T(0) || n <= degree + 1) {
        auto [knots, coeffs] = bsplineCoefficients(x, y, degree);
        // Residual computation
        T residual = T(0);
        for (std::size_t i = 0; i < n; ++i) {
            T diff = y[i] - bsplineEvaluate(
                std::span<const T>(knots), std::span<const T>(coeffs), degree, x[i]);
            residual += weights[i] * diff * diff;
        }
        return { std::move(knots), std::move(coeffs), degree, T(0), residual };
    }

    // Smoothing spline: adaptively increase the number of knots
    // Initial knot count: degree+1 boundary knots + a small number of interior knots
    std::size_t min_interior = 1;
    std::size_t max_interior = n - degree - 1;

    // Determine the knot count by binary search
    std::size_t best_interior = min_interior;
    std::vector<T> best_knots;
    std::vector<T> best_coeffs;
    T best_residual = std::numeric_limits<T>::max();

    auto try_fit = [&](std::size_t num_interior) -> T {
        // Build the knot vector
        std::size_t num_basis = num_interior + degree + 1;
        if (num_basis > n) num_basis = n;
        std::size_t num_knots = num_basis + degree + 1;
        std::vector<T> knots(num_knots);

        // Clamped endpoints
        for (std::size_t i = 0; i <= degree; ++i) {
            knots[i] = x.front();
            knots[num_knots - 1 - i] = x.back();
        }

        // Interior knots: place the data points at equal quantiles
        for (std::size_t j = 0; j < num_interior; ++j) {
            T frac = static_cast<T>(j + 1) / static_cast<T>(num_interior + 1);
            // Use the data quantile points
            T idx = frac * static_cast<T>(n - 1);
            std::size_t lo = static_cast<std::size_t>(idx);
            T t = idx - static_cast<T>(lo);
            if (lo >= n - 1) {
                knots[degree + 1 + j] = x[n - 1];
            } else {
                knots[degree + 1 + j] = x[lo] * (T(1) - t) + x[lo + 1] * t;
            }
        }

        // Least-squares fitting: collocation matrix A (n × num_basis)
        // Least-squares problem: min Σ w_i (y_i - Σ c_j B_j(x_i))²
        // -> normal equations: (A^T W A) c = A^T W y
        std::span<const T> knots_span(knots);

        Matrix<T> A(static_cast<typename Matrix<T>::size_type>(n),
                    static_cast<typename Matrix<T>::size_type>(num_basis));
        for (std::size_t i = 0; i < n; ++i) {
            for (std::size_t j = 0; j < num_basis; ++j) {
                A(static_cast<typename Matrix<T>::size_type>(i),
                  static_cast<typename Matrix<T>::size_type>(j)) =
                    bsplineBasis(degree, j, x[i], knots_span);
            }
        }

        // Build W^{1/2} A and W^{1/2} y
        Matrix<T> WA(static_cast<typename Matrix<T>::size_type>(n),
                     static_cast<typename Matrix<T>::size_type>(num_basis));
        Vector<T> Wy(static_cast<typename Vector<T>::size_type>(n));
        for (std::size_t i = 0; i < n; ++i) {
            T sw = std::sqrt(weights[i]);
            for (std::size_t j = 0; j < num_basis; ++j) {
                WA(static_cast<typename Matrix<T>::size_type>(i),
                   static_cast<typename Matrix<T>::size_type>(j)) =
                    sw * A(static_cast<typename Matrix<T>::size_type>(i),
                           static_cast<typename Matrix<T>::size_type>(j));
            }
            Wy[static_cast<typename Vector<T>::size_type>(i)] = sw * y[i];
        }

        // Normal equations: (WA)^T (WA) c = (WA)^T Wy
        Matrix<T> AtA(static_cast<typename Matrix<T>::size_type>(num_basis),
                      static_cast<typename Matrix<T>::size_type>(num_basis));
        Vector<T> Atb(static_cast<typename Vector<T>::size_type>(num_basis));
        for (std::size_t i = 0; i < num_basis; ++i) {
            for (std::size_t j = 0; j < num_basis; ++j) {
                T sum = T(0);
                for (std::size_t k = 0; k < n; ++k) {
                    sum += WA(static_cast<typename Matrix<T>::size_type>(k),
                              static_cast<typename Matrix<T>::size_type>(i))
                         * WA(static_cast<typename Matrix<T>::size_type>(k),
                              static_cast<typename Matrix<T>::size_type>(j));
                }
                AtA(static_cast<typename Matrix<T>::size_type>(i),
                    static_cast<typename Matrix<T>::size_type>(j)) = sum;
            }
            T sum = T(0);
            for (std::size_t k = 0; k < n; ++k) {
                sum += WA(static_cast<typename Matrix<T>::size_type>(k),
                          static_cast<typename Matrix<T>::size_type>(i))
                     * Wy[static_cast<typename Vector<T>::size_type>(k)];
            }
            Atb[static_cast<typename Vector<T>::size_type>(i)] = sum;
        }

        // Solve
        Vector<T> c = algorithms::solve(AtA, Atb, algorithms::SolverType::LU);

        // Residual computation
        T residual = T(0);
        std::vector<T> coeffs(num_basis);
        for (std::size_t i = 0; i < num_basis; ++i)
            coeffs[i] = c[static_cast<typename Vector<T>::size_type>(i)];
        for (std::size_t i = 0; i < n; ++i) {
            T val = bsplineEvaluate(knots_span, std::span<const T>(coeffs), degree, x[i]);
            T diff = y[i] - val;
            residual += weights[i] * diff * diff;
        }

        if (residual < best_residual || best_knots.empty()) {
            best_residual = residual;
            best_knots = knots;
            best_coeffs = std::move(coeffs);
            best_interior = num_interior;
        }

        return residual;
    };

    // Increase the number of knots, fitting until the residual ≤ s
    for (std::size_t ni = min_interior; ni <= max_interior; ++ni) {
        T res = try_fit(ni);
        if (res <= s) break;
        // Stop if the residual got worse than before (avoid overfitting)
        if (ni > min_interior + 2 && res > best_residual * T(1.1)) break;
    }

    return { std::move(best_knots), std::move(best_coeffs), degree, s, best_residual };
}

/// std::vector version of splineFit
template<concepts::OrderedField T>
[[nodiscard]]
SplineFitResult<T> splineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& w = {},
    T s = T(-1),
    std::size_t degree = 3)
{
    return splineFit(std::span<const T>(x), std::span<const T>(y),
                     w.empty() ? std::span<const T>{} : std::span<const T>(w),
                     s, degree);
}

/**
 * @brief Evaluate a fitted spline (equivalent to FITPACK splev)
 *
 * Evaluates at an arbitrary point using the result of splineFit.
 *
 * @param fit     Result of splineFit
 * @param x       Evaluation point
 * @return Spline value
 */
template<concepts::OrderedField T>
[[nodiscard]]
T splineEval(const SplineFitResult<T>& fit, T x) {
    return bsplineEvaluate(
        std::span<const T>(fit.knots),
        std::span<const T>(fit.coefficients),
        fit.degree, x);
}

/**
 * @brief Multi-point evaluation of a fitted spline
 */
template<concepts::OrderedField T>
[[nodiscard]]
std::vector<T> splineEval(const SplineFitResult<T>& fit, std::span<const T> xs) {
    std::vector<T> result(xs.size());
    for (std::size_t i = 0; i < xs.size(); ++i)
        result[i] = splineEval(fit, xs[i]);
    return result;
}

// =====================================================================
// Parametric spline
// =====================================================================

/**
 * @brief Result of parametric spline fitting
 */
template<concepts::OrderedField T>
struct ParametricSplineResult {
    SplineFitResult<T> x_fit;   ///< Spline fit of x(t)
    SplineFitResult<T> y_fit;   ///< Spline fit of y(t)
    std::vector<T> t;           ///< Parameter values
};

/**
 * @brief Parametric spline fitting
 *
 * Spline-fits a planar curve (x(t), y(t)) using the arc-length parameter.
 *
 * @param x       x-coordinates of the data points
 * @param y       y-coordinates of the data points
 * @param s       Smoothing parameter (0 = interpolation, < 0 = automatic)
 * @param degree  Spline degree (default 3)
 * @return ParametricSplineResult<T>
 */
template<concepts::OrderedField T>
[[nodiscard]]
ParametricSplineResult<T> parametricSplineFit(
    std::span<const T> x,
    std::span<const T> y,
    T s = T(-1),
    std::size_t degree = 3)
{
    std::size_t n = x.size();
    if (n < 2)
        throw std::invalid_argument("parametricSplineFit: need at least 2 points");
    if (y.size() != n)
        throw std::invalid_argument("parametricSplineFit: x and y must have same size");

    // Compute the arc-length parameter t (cumulative chord length)
    std::vector<T> t(n);
    t[0] = T(0);
    for (std::size_t i = 1; i < n; ++i) {
        T dx = x[i] - x[i - 1];
        T dy = y[i] - y[i - 1];
        t[i] = t[i - 1] + std::sqrt(dx * dx + dy * dy);
    }
    // Normalize to [0, 1]
    if (t.back() > T(0)) {
        T total = t.back();
        for (auto& ti : t) ti /= total;
    }

    // Fit x(t) and y(t) separately
    std::span<const T> t_span(t);
    auto x_fit = splineFit(t_span, x, std::span<const T>{}, s, degree);
    auto y_fit = splineFit(t_span, y, std::span<const T>{}, s, degree);

    return { std::move(x_fit), std::move(y_fit), std::move(t) };
}

/// std::vector version of parametricSplineFit
template<concepts::OrderedField T>
[[nodiscard]]
ParametricSplineResult<T> parametricSplineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T s = T(-1),
    std::size_t degree = 3)
{
    return parametricSplineFit(std::span<const T>(x), std::span<const T>(y), s, degree);
}

/**
 * @brief Evaluate a parametric spline
 */
template<concepts::OrderedField T>
[[nodiscard]]
std::pair<T, T> parametricSplineEval(const ParametricSplineResult<T>& fit, T t) {
    return { splineEval(fit.x_fit, t), splineEval(fit.y_fit, t) };
}

// =====================================================================
// 2D tensor-product spline (surface fitting)
// =====================================================================

/**
 * @brief Result of 2D tensor-product spline fitting
 */
template<concepts::OrderedField T>
struct SurfaceSplineResult {
    std::vector<T> knots_x;     ///< Knot vector in the x direction
    std::vector<T> knots_y;     ///< Knot vector in the y direction
    std::vector<T> coefficients;///< Coefficients (nx × ny flat array, row-major)
    std::size_t nx;             ///< Number of basis functions in the x direction
    std::size_t ny;             ///< Number of basis functions in the y direction
    std::size_t degree;         ///< Spline degree
};

/**
 * @brief 2D tensor-product spline fitting
 *
 * Tensor-product B-spline fitting of grid data z[i][j] = f(x[i], y[j]).
 *
 * S(x,y) = ΣΣ c_{ij} B_i(x) B_j(y)
 *
 * @param x       Grid coordinates in the x direction (mx of them, ascending)
 * @param y       Grid coordinates in the y direction (my of them, ascending)
 * @param z       Value matrix z[i*my + j] = f(x[i], y[j]) (mx*my, row-major)
 * @param degree  Spline degree (default 3)
 * @return SurfaceSplineResult<T>
 */
template<concepts::OrderedField T>
[[nodiscard]]
SurfaceSplineResult<T> surfaceSplineFit(
    std::span<const T> x,
    std::span<const T> y,
    std::span<const T> z,
    std::size_t degree = 3)
{
    std::size_t mx = x.size();
    std::size_t my = y.size();
    if (mx < degree + 1 || my < degree + 1)
        throw std::invalid_argument("surfaceSplineFit: need at least (degree+1) points in each direction");
    if (z.size() != mx * my)
        throw std::invalid_argument("surfaceSplineFit: z must have mx*my elements");

    // Knot vector in the x direction (clamped)
    std::size_t nkx = mx + degree + 1;
    std::vector<T> knots_x(nkx);
    for (std::size_t i = 0; i <= degree; ++i) {
        knots_x[i] = x.front();
        knots_x[nkx - 1 - i] = x.back();
    }
    for (std::size_t j = 1; j < mx - degree; ++j) {
        T sum = T(0);
        for (std::size_t i = j; i < j + degree; ++i) sum += x[i];
        knots_x[j + degree] = sum / static_cast<T>(degree);
    }

    // Knot vector in the y direction
    std::size_t nky = my + degree + 1;
    std::vector<T> knots_y(nky);
    for (std::size_t i = 0; i <= degree; ++i) {
        knots_y[i] = y.front();
        knots_y[nky - 1 - i] = y.back();
    }
    for (std::size_t j = 1; j < my - degree; ++j) {
        T sum = T(0);
        for (std::size_t i = j; i < j + degree; ++i) sum += y[i];
        knots_y[j + degree] = sum / static_cast<T>(degree);
    }

    // Tensor-product decomposition: for each x row, find the y-direction B-spline coefficients,
    // then for each y column of those coefficients, find the x-direction B-spline coefficients

    // Step 1: for each x row i, B-spline fit z[i,:] in the y direction
    // -> alpha[i][j] (mx × my)
    std::span<const T> ky_span(knots_y);

    // y-direction collocation matrix By (my × my)
    Matrix<T> By(static_cast<typename Matrix<T>::size_type>(my),
                 static_cast<typename Matrix<T>::size_type>(my));
    for (std::size_t i = 0; i < my; ++i) {
        for (std::size_t j = 0; j < my; ++j) {
            By(static_cast<typename Matrix<T>::size_type>(i),
               static_cast<typename Matrix<T>::size_type>(j)) =
                bsplineBasis(degree, j, y[i], ky_span);
        }
    }
    // Boundary rows
    for (std::size_t j = 0; j < my; ++j) {
        By(0, static_cast<typename Matrix<T>::size_type>(j)) = T(0);
        By(static_cast<typename Matrix<T>::size_type>(my - 1),
           static_cast<typename Matrix<T>::size_type>(j)) = T(0);
    }
    By(0, 0) = T(1);
    By(static_cast<typename Matrix<T>::size_type>(my - 1),
       static_cast<typename Matrix<T>::size_type>(my - 1)) = T(1);

    // alpha[i][j]: y-direction coefficients for mx rows
    std::vector<T> alpha(mx * my);
    for (std::size_t i = 0; i < mx; ++i) {
        Vector<T> rhs(static_cast<typename Vector<T>::size_type>(my));
        for (std::size_t j = 0; j < my; ++j)
            rhs[static_cast<typename Vector<T>::size_type>(j)] = z[i * my + j];
        Vector<T> c = algorithms::solve(By, rhs, algorithms::SolverType::LU);
        for (std::size_t j = 0; j < my; ++j)
            alpha[i * my + j] = c[static_cast<typename Vector<T>::size_type>(j)];
    }

    // Step 2: for each y column j, B-spline fit alpha[:,j] in the x direction
    // -> coefficients[i*my + j] (mx × my)
    std::span<const T> kx_span(knots_x);

    Matrix<T> Bx(static_cast<typename Matrix<T>::size_type>(mx),
                 static_cast<typename Matrix<T>::size_type>(mx));
    for (std::size_t i = 0; i < mx; ++i) {
        for (std::size_t j = 0; j < mx; ++j) {
            Bx(static_cast<typename Matrix<T>::size_type>(i),
               static_cast<typename Matrix<T>::size_type>(j)) =
                bsplineBasis(degree, j, x[i], kx_span);
        }
    }
    for (std::size_t j = 0; j < mx; ++j) {
        Bx(0, static_cast<typename Matrix<T>::size_type>(j)) = T(0);
        Bx(static_cast<typename Matrix<T>::size_type>(mx - 1),
           static_cast<typename Matrix<T>::size_type>(j)) = T(0);
    }
    Bx(0, 0) = T(1);
    Bx(static_cast<typename Matrix<T>::size_type>(mx - 1),
       static_cast<typename Matrix<T>::size_type>(mx - 1)) = T(1);

    std::vector<T> coefficients(mx * my);
    for (std::size_t j = 0; j < my; ++j) {
        Vector<T> rhs(static_cast<typename Vector<T>::size_type>(mx));
        for (std::size_t i = 0; i < mx; ++i)
            rhs[static_cast<typename Vector<T>::size_type>(i)] = alpha[i * my + j];
        Vector<T> c = algorithms::solve(Bx, rhs, algorithms::SolverType::LU);
        for (std::size_t i = 0; i < mx; ++i)
            coefficients[i * my + j] = c[static_cast<typename Vector<T>::size_type>(i)];
    }

    return { std::move(knots_x), std::move(knots_y),
             std::move(coefficients), mx, my, degree };
}

/// std::vector version of surfaceSplineFit
template<concepts::OrderedField T>
[[nodiscard]]
SurfaceSplineResult<T> surfaceSplineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& z,
    std::size_t degree = 3)
{
    return surfaceSplineFit(std::span<const T>(x), std::span<const T>(y),
                            std::span<const T>(z), degree);
}

/**
 * @brief Evaluate a 2D tensor-product spline
 */
template<concepts::OrderedField T>
[[nodiscard]]
T surfaceSplineEval(const SurfaceSplineResult<T>& fit, T x, T y) {
    std::span<const T> kx(fit.knots_x);
    std::span<const T> ky(fit.knots_y);
    T result = T(0);
    for (std::size_t i = 0; i < fit.nx; ++i) {
        T bx = bsplineBasis(fit.degree, i, x, kx);
        if (bx == T(0)) continue;
        for (std::size_t j = 0; j < fit.ny; ++j) {
            T by = bsplineBasis(fit.degree, j, y, ky);
            result += fit.coefficients[i * fit.ny + j] * bx * by;
        }
    }
    return result;
}

// =====================================================================
// Nearest-neighbor interpolation — 2D
// =====================================================================

/**
 * @brief 2D nearest-neighbor interpolation
 *
 * Returns the value of the nearest grid point on the grid as-is.
 * Used for image downscaling and pixel-art upscaling.
 *
 * @param grid  Row-major 2D data (rows × cols)
 * @param rows  Number of rows in the grid
 * @param cols  Number of columns in the grid
 * @param x     x-coordinate of the evaluation point (0-based, column direction)
 * @param y     y-coordinate of the evaluation point (0-based, row direction)
 * @return Value at the nearest grid point
 */
template<concepts::OrderedField T>
[[nodiscard]] T nearestNeighbor2D(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y) {

    // Nearest-neighbor index (rounding)
    auto ix = static_cast<std::size_t>(std::max(T(0), std::min(
        std::round(x), static_cast<T>(cols - 1))));
    auto iy = static_cast<std::size_t>(std::max(T(0), std::min(
        std::round(y), static_cast<T>(rows - 1))));

    return grid[iy * cols + ix];
}

// =====================================================================
// Bicubic interpolation — 2D
// =====================================================================

/**
 * @brief Bicubic interpolation
 *
 * Uses 4×4 = 16 neighboring grid points to interpolate smoothly with a cubic polynomial.
 * Equivalent to OpenCV INTER_CUBIC and Photoshop's bicubic method.
 * Uses the Catmull-Rom spline kernel (a = -0.5).
 *
 * @param grid  Row-major 2D data (rows × cols)
 * @param rows  Number of rows in the grid
 * @param cols  Number of columns in the grid
 * @param x     x-coordinate of the evaluation point (0-based, column direction)
 * @param y     y-coordinate of the evaluation point (0-based, row direction)
 * @return Interpolated value
 */
template<concepts::OrderedField T>
[[nodiscard]] T bicubicInterpolate(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y) {

    // Catmull-Rom kernel (a = -0.5)
    auto kernel = [](T t) -> T {
        T at = std::abs(t);
        if (at <= T(1)) {
            return (T(3) / T(2)) * at * at * at
                 - (T(5) / T(2)) * at * at + T(1);
        } else if (at < T(2)) {
            return -(T(1) / T(2)) * at * at * at
                 + (T(5) / T(2)) * at * at
                 - T(4) * at + T(2);
        }
        return T(0);
    };

    // Clamped index
    auto clamp_idx = [](int v, int max_v) -> std::size_t {
        return static_cast<std::size_t>(std::max(0, std::min(v, max_v)));
    };

    int ix = static_cast<int>(detail::to_int64(std::floor(x)));
    int iy = static_cast<int>(detail::to_int64(std::floor(y)));
    T fx = x - static_cast<T>(ix);
    T fy = y - static_cast<T>(iy);

    int max_col = static_cast<int>(cols) - 1;
    int max_row = static_cast<int>(rows) - 1;

    T result = T(0);
    for (int m = -1; m <= 2; ++m) {
        T wy = kernel(fy - static_cast<T>(m));
        std::size_t row = clamp_idx(iy + m, max_row);
        for (int n = -1; n <= 2; ++n) {
            T wx = kernel(fx - static_cast<T>(n));
            std::size_t col = clamp_idx(ix + n, max_col);
            result += wy * wx * grid[row * cols + col];
        }
    }
    return result;
}

// =====================================================================
// Lanczos interpolation — 2D
// =====================================================================

/**
 * @brief 2D Lanczos interpolation
 *
 * High-quality sinc-based resampling. The Lanczos-a kernel:
 *   L(x) = sinc(x) * sinc(x/a)   when |x| < a
 *   L(x) = 0                     otherwise
 *
 * a = 2 (Lanczos-2, 4×4 window) or a = 3 (Lanczos-3, 6×6 window).
 * The default is a = 3 (highest quality).
 *
 * @param grid  Row-major 2D data (rows × cols)
 * @param rows  Number of rows in the grid
 * @param cols  Number of columns in the grid
 * @param x     x-coordinate of the evaluation point (0-based, column direction)
 * @param y     y-coordinate of the evaluation point (0-based, row direction)
 * @param a     Lanczos parameter (window width, default 3)
 * @return Interpolated value
 */
template<concepts::OrderedField T>
[[nodiscard]] T lanczosInterpolate2D(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y,
    int a = 3) {

    // Lanczos kernel: compute the weights (which directly affect the interpolated result) in T.
    // Going through double would lose the precision of Float (use the generic_math helpers).
    const T pi = detail::generic_pi<T>();
    const T da = T(a);
    auto lanczos_kernel = [&](T t) -> T {
        if (t == T(0)) return T(1);
        if (detail::generic_abs(t) >= da) return T(0);
        T pi_t = pi * t;
        return da * detail::generic_sin(pi_t) * detail::generic_sin(pi_t / da)
                  / (pi_t * pi_t);
    };

    auto clamp_idx = [](int v, int max_v) -> std::size_t {
        return static_cast<std::size_t>(std::max(0, std::min(v, max_v)));
    };

    int ix = static_cast<int>(detail::to_int64(std::floor(x)));
    int iy = static_cast<int>(detail::to_int64(std::floor(y)));
    T fx = x - static_cast<T>(ix);
    T fy = y - static_cast<T>(iy);

    int max_col = static_cast<int>(cols) - 1;
    int max_row = static_cast<int>(rows) - 1;

    T result = T(0);
    T weight_sum = T(0);

    for (int m = -(a - 1); m <= a; ++m) {
        T wy = lanczos_kernel(fy - static_cast<T>(m));
        std::size_t row = clamp_idx(iy + m, max_row);
        for (int n = -(a - 1); n <= a; ++n) {
            T wx = lanczos_kernel(fx - static_cast<T>(n));
            std::size_t col = clamp_idx(ix + n, max_col);
            T w = wy * wx;
            result += w * grid[row * cols + col];
            weight_sum += w;
        }
    }

    // Weight normalization (energy conservation near the boundary)
    if (weight_sum != T(0)) {
        result /= weight_sum;
    }
    return result;
}

} // namespace sangi

#endif // INTERPOLATION_HPP