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

// root_finding_1d.hpp
#ifndef SANGI_ROOT_FINDING_1D_HPP
#define SANGI_ROOT_FINDING_1D_HPP

#include "root_finding_base.hpp"
#include <math/core/Complex.hpp>
#include <vector>

namespace sangi {

/**
 * @brief Root finding by the bisection method
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param a Lower bound of the interval
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> bisection(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);

    // Check whether a root lies in the interval [a, b]
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(0, std::abs(b - a));
    }

    // If fa == 0 or fb == 0, that point is the root
    if (criteria.is_function_converged(fa))
        return RootFindingResult<T>::success(a, 0, std::numeric_limits<T>::epsilon());

    if (criteria.is_function_converged(fb))
        return RootFindingResult<T>::success(b, 0, std::numeric_limits<T>::epsilon());

    T c_prev = a;
    size_t iterations = 0;

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        T c = (a + b) / 2;  // Compute the midpoint
        T fc = f(c);

        // Combined convergence check
        if (criteria.is_function_converged(fc) || criteria.is_interval_converged(a, b)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(b - a) / 2);
        }

        // Narrow the interval down to half
        if (fa * fc < 0) {
            b = c;
            fb = fc;
        }
        else {
            a = c;
            fa = fc;
        }

        // Treat as converged also when the change from the previous midpoint is negligible
        if (criteria.is_variable_converged(c, c_prev)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(c - c_prev));
        }

        c_prev = c;
    }

    // Maximum number of iterations reached; return the best approximation
    T c = (a + b) / 2;
    return RootFindingResult<T>::partial_success(c, false, iterations, std::abs(b - a) / 2);
}

/**
 * @brief Root finding by the Newton-Raphson method
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param df Derivative of the function f
 * @param x0 Initial value
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> newton_raphson(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    T x_prev = x0;
    T fx_prev = std::numeric_limits<T>::max();
    size_t iterations = 0;

    for (iterations = 0; iterations < criteria.max_iterations; ++iterations) {
        T fx = f(x);

        // Convergence check by function value
        if (criteria.is_function_converged(fx, fx_prev)) {
            return RootFindingResult<T>::success(x, iterations, std::abs(x - x_prev));
        }

        T dfx = df(x);

        // When the derivative is close to zero
        if (std::abs(dfx) < std::numeric_limits<T>::epsilon() * 10) {
            // Treat as converged if the current function value is small enough
            if (criteria.is_function_converged(fx)) {
                return RootFindingResult<T>::success(x, iterations, std::abs(fx) / (std::abs(dfx) + std::numeric_limits<T>::epsilon()));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx));
        }

        x_prev = x;
        T x_new = x - fx / dfx;

        // Convergence check by variable value
        if (criteria.is_variable_converged(x_new, x_prev)) {
            return RootFindingResult<T>::success(x_new, iterations + 1, std::abs(x_new - x_prev));
        }

        x = x_new;
        fx_prev = fx;
    }

    // When the maximum number of iterations is reached
    return RootFindingResult<T>::failure(iterations, std::abs(f(x)));
}

/**
 * @brief Root finding by the secant method
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param x0 First initial value
 * @param x1 Second initial value
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> secant_method(
    const std::function<T(T)>& f,
    T x0, T x1,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fx0 = f(x0);
    size_t iterations = 0;

    // When x0 is already a root
    if (criteria.is_function_converged(fx0)) {
        return RootFindingResult<T>::success(x0, iterations, std::numeric_limits<T>::epsilon());
    }

    T fx1 = f(x1);

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        // Convergence check by function value
        if (criteria.is_function_converged(fx1, fx0)) {
            return RootFindingResult<T>::success(x1, iterations, std::abs(x1 - x0));
        }

        // When the slope of the secant is close to zero
        if (std::abs(fx1 - fx0) < std::numeric_limits<T>::epsilon()) {
            // If the slope is nearly zero, treat as converged when fx1 is small enough
            if (criteria.is_function_converged(fx1)) {
                return RootFindingResult<T>::success(x1, iterations, std::abs(fx1));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx1));
        }

        // Update formula of the secant method
        T x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0);

        // Convergence check by variable value
        if (criteria.is_variable_converged(x2, x1)) {
            return RootFindingResult<T>::success(x2, iterations, std::abs(x2 - x1));
        }

        // Update the values
        x0 = x1;
        fx0 = fx1;
        x1 = x2;
        fx1 = f(x1);
    }

    // Maximum number of iterations reached without convergence
    return RootFindingResult<T>::failure(iterations, std::abs(fx1));
}

/**
 * @brief Root finding by Brent's method
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param a Lower bound of the interval
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> brent_method(
    const std::function<T(T)>& f,
    T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);
    size_t iterations = 0;

    // Check whether a root lies in the interval [a, b]
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(iterations, std::abs(b - a));
    }

    // When fa == 0 or fb == 0
    if (std::abs(fa) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(a, iterations, std::numeric_limits<T>::epsilon());
    }
    if (std::abs(fb) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(b, iterations, std::numeric_limits<T>::epsilon());
    }

    // Swap a and b so that |fa| > |fb|
    if (std::abs(fa) < std::abs(fb)) {
        std::swap(a, b);
        std::swap(fa, fb);
    }

    T c = a;
    T fc = fa;
    T d = c;
    bool mflag = true;
    T s = b;  // Start the initial approximation from b
    T fs = fb;
    T s_prev = a;  // Set the previous value to a clearly different value

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        // Convergence check
        if (std::abs(fs) < criteria.abs_ftol && iterations > 1) {
            return RootFindingResult<T>::success(s, iterations, std::abs(s - s_prev));
        }

        // Convergence check by interval width
        if (criteria.is_interval_converged(a, b) && iterations > 1) {
            return RootFindingResult<T>::success(s, iterations, std::abs(b - a));
        }

        // Convergence check by variable value
        if (iterations > 1 && criteria.is_variable_converged(s, s_prev)) {
            return RootFindingResult<T>::success(s, iterations, std::abs(s - s_prev));
        }

        // Compute the next approximation point
        s_prev = s;
        T s_new;

        // Inverse quadratic interpolation or secant method calculation
        if (std::abs(fc - fa) > criteria.abs_ftol &&
            std::abs(fb - fc) > criteria.abs_ftol) {
            // Inverse quadratic interpolation
            s_new = a * fb * fc / ((fa - fb) * (fa - fc)) +
                b * fa * fc / ((fb - fa) * (fb - fc)) +
                c * fa * fb / ((fc - fa) * (fc - fb));
        }
        else {
            // Secant method
            s_new = b - fb * (b - a) / (fb - fa);
        }

        // Use bisection if s_new is out of bounds or not acceptable
        bool use_bisection =
            (s_new <= (3 * a + b) / 4 || s_new >= b) ||
            (mflag && std::abs(s_new - b) >= std::abs(b - c) / 2) ||
            (!mflag && std::abs(s_new - b) >= std::abs(c - d) / 2) ||
            (mflag && std::abs(b - c) < criteria.abs_xtol) ||
            (!mflag && std::abs(c - d) < criteria.abs_xtol);

        if (use_bisection) {
            s_new = (a + b) / 2;  // Bisection
            mflag = true;
        }
        else {
            mflag = false;
        }

        // Function evaluation at the new approximation
        s = s_new;
        fs = f(s);

        // Update the points
        d = c;
        c = b;
        fc = fb;

        // Update to the new interval that contains the root
        if (fa * fs < 0) {
            b = s;
            fb = fs;
        }
        else {
            a = s;
            fa = fs;
        }

        // Swap the points so that |fa| >= |fb|
        if (std::abs(fa) < std::abs(fb)) {
            std::swap(a, b);
            std::swap(fa, fb);
        }
    }

    // When the maximum number of iterations is reached
    // Return the point with the smallest function value
    if (std::abs(fa) < std::abs(fb)) {
        s = a;
        fs = fa;
    }
    else {
        s = b;
        fs = fb;
    }

    bool converged = std::abs(fs) < criteria.abs_ftol;
    return RootFindingResult<T>::partial_success(s, converged, iterations, std::abs(fs));
}

/**
 * @brief Root finding by the ITP method (Interpolate–Truncate–Project)
 *
 * The bracketing method of Oliveira & Takahashi (2020) "An Enhancement of the
 * Bisection Method Average Performance Preserving Minmax Optimality" (ACM TOMS
 * 47(1)). Each iteration determines the next point in three stages: (1) secant
 * interpolation x_f, (2) truncation that perturbs the interpolation point toward
 * the midpoint by the minimal amount δ=κ₁(b−a)^κ₂, and (3) projection toward the
 * center of the interval so as not to exceed the worst-case bisection iteration
 * count +n₀. It combines super-linear convergence speed with minmax optimality
 * against bisection (at worst bisection +n₀ iterations).
 *
 * @param f The function whose root is sought (requires f(a)·f(b) < 0)
 * @param a Lower bound of the interval
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria (uses abs_xtol as the half-interval tolerance ε)
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> itp_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);

    if (std::abs(fa) < criteria.abs_ftol)
        return RootFindingResult<T>::success(a, 0, std::numeric_limits<T>::epsilon());
    if (std::abs(fb) < criteria.abs_ftol)
        return RootFindingResult<T>::success(b, 0, std::numeric_limits<T>::epsilon());
    if (fa * fb > T(0))
        return RootFindingResult<T>::failure(0, std::abs(b - a));

    // Half-interval tolerance ε. Reuses abs_xtol (invalid values fall back to a machine-ε basis).
    T eps = criteria.abs_xtol;
    if (eps <= T(0)) eps = std::numeric_limits<T>::epsilon() * 100;

    const T width0 = std::abs(b - a);
    // Scale-invariant default parameters. κ₂=2 satisfying κ₂ ∈ [1, 1+φ), and κ₁ normalized by the initial interval width.
    const T k1 = T(0.2) / width0;
    const T k2 = T(2);
    const int n0 = 1;

    int n_half = static_cast<int>(std::ceil(std::log2(width0 / (T(2) * eps))));
    if (n_half < 0) n_half = 0;
    const int n_max = n_half + n0;

    size_t iterations = 0;
    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        const T x_half = (a + b) / T(2);

        if (std::abs(b - a) / T(2) <= eps)
            return RootFindingResult<T>::success(x_half, iterations, std::abs(b - a) / T(2));

        // (1) Interpolation: secant (regula falsi) estimate point
        const T denom = fa - fb;
        const T x_f = (std::abs(denom) > std::numeric_limits<T>::min())
                          ? (b * fa - a * fb) / denom
                          : x_half;

        // (2) Truncation: truncate x_f toward the midpoint by the minimal perturbation δ
        const T sigma = (x_half - x_f >= T(0)) ? T(1) : T(-1);
        const T delta = k1 * std::pow(std::abs(b - a), k2);
        const T x_t = (delta <= std::abs(x_half - x_f)) ? (x_f + sigma * delta) : x_half;

        // (3) Projection: project to within distance r of the midpoint (guarantees the worst-case bisection iteration count)
        const int jexp = n_max - static_cast<int>(iterations - 1);
        T r = eps * std::pow(T(2), T(jexp)) - std::abs(b - a) / T(2);
        if (r < T(0)) r = T(0);
        const T x_itp = (std::abs(x_t - x_half) <= r) ? x_t : (x_half - sigma * r);

        const T y = f(x_itp);
        if (std::abs(y) < criteria.abs_ftol)
            return RootFindingResult<T>::success(x_itp, iterations, std::abs(b - a) / T(2));

        // Interval update (fa, fb always keep opposite signs)
        if ((y > T(0)) == (fb > T(0))) { b = x_itp; fb = y; }
        else                          { a = x_itp; fa = y; }
    }

    const T x_final = (a + b) / T(2);
    return RootFindingResult<T>::partial_success(x_final, false, iterations, std::abs(b - a) / T(2));
}

/**
 * @brief Root finding by Ridders' method
 *
 * A kind of bracketing method. It applies an exponential transformation to the
 * bisection midpoint to achieve quadratic convergence.
 * x_new = x_mid + (x_mid - a) * sign(fa - fb) * f_mid / sqrt(f_mid² - fa*fb)
 *
 * @param f The function whose root is sought
 * @param a Lower bound of the interval (requires f(a)·f(b) < 0)
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> ridders_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);
    size_t iterations = 0;

    // Check whether a root lies in the interval [a, b]
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(iterations, std::abs(b - a));
    }

    if (std::abs(fa) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(a, iterations, std::numeric_limits<T>::epsilon());
    }
    if (std::abs(fb) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(b, iterations, std::numeric_limits<T>::epsilon());
    }

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        T mid = (a + b) / 2;
        T fm = f(mid);

        // Ridders' formula
        T disc = fm * fm - fa * fb;
        if (disc < T(0)) disc = T(0);
        T sq = std::sqrt(disc);
        if (sq == T(0)) {
            return RootFindingResult<T>::success(mid, iterations, std::abs(b - a) / 2);
        }

        T sign = (fa - fb > T(0)) ? T(1) : T(-1);
        T x_new = mid + (mid - a) * sign * fm / sq;

        T fx_new = f(x_new);

        // Convergence check
        if (std::abs(fx_new) < criteria.abs_ftol) {
            return RootFindingResult<T>::success(x_new, iterations, std::abs(fx_new));
        }

        // Bracket update
        if (fm * fx_new < T(0)) {
            // A root lies between mid and x_new
            if (mid < x_new) {
                a = mid; fa = fm;
                b = x_new; fb = fx_new;
            } else {
                a = x_new; fa = fx_new;
                b = mid; fb = fm;
            }
        } else if (fa * fx_new < T(0)) {
            b = x_new; fb = fx_new;
        } else {
            a = x_new; fa = fx_new;
        }

        // Convergence check by interval width
        if (criteria.is_interval_converged(a, b)) {
            return RootFindingResult<T>::success((a + b) / 2, iterations, std::abs(b - a));
        }
    }

    T best = (std::abs(fa) < std::abs(fb)) ? a : b;
    T best_f = std::min(std::abs(fa), std::abs(fb));
    return RootFindingResult<T>::partial_success(best, false, iterations, best_f);
}

/**
 * @brief Root finding by Regula Falsi (false position method, basic form)
 *
 * The next trial value is the intersection of the x-axis with the line joining
 * the two ends of the interval [a, b] (f(a)·f(b) < 0); the endpoint with the same
 * sign as f(c) is replaced by c, keeping a sign-change bracketing interval. Since
 * the interval always brackets a root, convergence is guaranteed, but when the
 * function is convex (or concave) one endpoint may stay fixed, causing stagnation
 * that can degrade convergence to linear.
 *
 * This function is the basic form without stagnation avoidance (for reference and
 * educational use). For practical use, the improved versions that avoid stagnation
 * are recommended: illinois_method (Illinois modification), king_method (Pegasus
 * family), and anderson_bjork_method.
 *
 * @param f The function whose root is sought
 * @param a Lower bound of the interval (requires f(a)·f(b) < 0)
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> regula_falsi(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);
    size_t iterations = 0;

    // Check whether a root lies in the interval [a, b]
    if (fa * fb > T(0)) {
        return RootFindingResult<T>::failure(iterations, std::abs(b - a));
    }
    if (std::abs(fa) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(a, iterations, std::numeric_limits<T>::epsilon());
    }
    if (std::abs(fb) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(b, iterations, std::numeric_limits<T>::epsilon());
    }

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        // Linear interpolation
        T c = (a * fb - b * fa) / (fb - fa);
        T fc = f(c);

        // Convergence check
        if (std::abs(fc) < criteria.abs_ftol) {
            return RootFindingResult<T>::success(c, iterations, std::abs(fc));
        }
        if (criteria.is_interval_converged(a, b)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(b - a));
        }

        // Bracket update (no stagnation avoidance = plain Regula Falsi)
        if (fa * fc < T(0)) {
            b = c; fb = fc;
        } else {
            a = c; fa = fc;
        }
    }

    T best = (std::abs(fa) < std::abs(fb)) ? a : b;
    T best_f = std::min(std::abs(fa), std::abs(fb));
    return RootFindingResult<T>::partial_success(best, false, iterations, best_f);
}

/**
 * @brief Root finding by the Illinois method (stagnation-avoiding false position)
 *
 * An improved version of Regula Falsi (false position method). It computes the next
 * approximation by linear interpolation, and on same-side stagnation applies the
 * Illinois modification (halving the weight) to avoid stagnation and guarantee
 * convergence. For the basic form without stagnation avoidance, see regula_falsi.
 *
 * @param f The function whose root is sought
 * @param a Lower bound of the interval (requires f(a)·f(b) < 0)
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> illinois_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);
    size_t iterations = 0;

    // Check whether a root lies in the interval [a, b]
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(iterations, std::abs(b - a));
    }

    if (std::abs(fa) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(a, iterations, std::numeric_limits<T>::epsilon());
    }
    if (std::abs(fb) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(b, iterations, std::numeric_limits<T>::epsilon());
    }

    int side = 0;  // 0=initial, 1=a-side changed, -1=b-side changed (for the Illinois modification)

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        // Linear interpolation
        T c = (a * fb - b * fa) / (fb - fa);
        T fc = f(c);

        // Convergence check
        if (std::abs(fc) < criteria.abs_ftol) {
            return RootFindingResult<T>::success(c, iterations, std::abs(fc));
        }

        if (criteria.is_interval_converged(a, b)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(b - a));
        }

        // Bracket update + Illinois modification
        if (fa * fc < T(0)) {
            b = c; fb = fc;
            if (side == 1) {
                // The a-side was not changed twice in a row → Illinois modification
                fa /= T(2);
            }
            side = 1;
        } else {
            a = c; fa = fc;
            if (side == -1) {
                fb /= T(2);
            }
            side = -1;
        }
    }

    T best = (std::abs(fa) < std::abs(fb)) ? a : b;
    T best_f = std::min(std::abs(fa), std::abs(fb));
    return RootFindingResult<T>::partial_success(best, false, iterations, best_f);
}

/**
 * @brief Root finding by Halley's method (cubic convergence)
 *
 * Uses f, f', f'' to achieve higher-order convergence than Newton's method.
 * Update formula: x_{n+1} = x_n - 2·f·f' / (2·f'² - f·f'')
 *
 * @param f   The function whose root is sought
 * @param df  First derivative of f
 * @param d2f Second derivative of f
 * @param x0  Initial value
 * @param criteria Convergence criteria
 */
template<concepts::OrderedField T>
RootFindingResult<T> halley_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    T x_prev = x0;
    size_t iterations = 0;

    for (iterations = 0; iterations < criteria.max_iterations; ++iterations) {
        T fx = f(x);

        if (criteria.is_function_converged(fx)) {
            return RootFindingResult<T>::success(x, iterations, std::abs(x - x_prev));
        }

        T dfx = df(x);
        T d2fx = d2f(x);

        // Denominator: 2·f'² - f·f''
        T denom = T(2) * dfx * dfx - fx * d2fx;
        if (std::abs(denom) < std::numeric_limits<T>::epsilon() * T(10)) {
            if (criteria.is_function_converged(fx)) {
                return RootFindingResult<T>::success(x, iterations, std::abs(fx));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx));
        }

        x_prev = x;
        T x_new = x - T(2) * fx * dfx / denom;

        if (criteria.is_variable_converged(x_new, x_prev)) {
            return RootFindingResult<T>::success(x_new, iterations + 1, std::abs(x_new - x_prev));
        }

        x = x_new;
    }

    return RootFindingResult<T>::failure(iterations, std::abs(f(x)));
}

/**
 * @brief Root finding by Schröder's method (cubic convergence, Householder order 2)
 *
 * The order-2 form of Householder's method using f, f', f''.
 * It maintains quadratic convergence even for multiple roots.
 * Update formula: x_{n+1} = x_n - f·f' / (f'² - f·f'')
 *
 * @param f   The function whose root is sought
 * @param df  First derivative of f
 * @param d2f Second derivative of f
 * @param x0  Initial value
 * @param criteria Convergence criteria
 */
template<concepts::OrderedField T>
RootFindingResult<T> schroder_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    T x_prev = x0;
    size_t iterations = 0;

    for (iterations = 0; iterations < criteria.max_iterations; ++iterations) {
        T fx = f(x);

        if (criteria.is_function_converged(fx)) {
            return RootFindingResult<T>::success(x, iterations, std::abs(x - x_prev));
        }

        T dfx = df(x);
        T d2fx = d2f(x);

        // Denominator: f'² - f·f''
        T denom = dfx * dfx - fx * d2fx;
        if (std::abs(denom) < std::numeric_limits<T>::epsilon() * T(10)) {
            if (criteria.is_function_converged(fx)) {
                return RootFindingResult<T>::success(x, iterations, std::abs(fx));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx));
        }

        x_prev = x;
        T x_new = x - fx * dfx / denom;

        if (criteria.is_variable_converged(x_new, x_prev)) {
            return RootFindingResult<T>::success(x_new, iterations + 1, std::abs(x_new - x_prev));
        }

        x = x_new;
    }

    return RootFindingResult<T>::failure(iterations, std::abs(f(x)));
}

/**
 * @brief Root finding by the Anderson-Björk method
 *
 * An improved version of Regula Falsi (false position method) that converges faster
 * than Illinois (weight 1/2). On same-side stagnation it multiplies fa by the weight
 * m = 1 - fc/fb (m = 1/2 if m <= 0) to accelerate convergence.
 *
 * Reference: Anderson & Björck, "A new high order method of regula falsi type for
 *       computing a root of an equation", BIT (1973).
 *
 * @param f The function whose root is sought
 * @param a Lower bound of the interval (requires f(a)·f(b) < 0)
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> anderson_bjork_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);
    size_t iterations = 0;

    // Check whether a root lies in the interval [a, b]
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(0, std::abs(b - a));
    }

    if (std::abs(fa) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(a, 0, std::numeric_limits<T>::epsilon());
    }
    if (std::abs(fb) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(b, 0, std::numeric_limits<T>::epsilon());
    }

    int side = 0;  // 0=initial, 1=a-side changed, -1=b-side changed (for the Anderson-Björk modification)

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        // Linear interpolation
        T c = (a * fb - b * fa) / (fb - fa);
        T fc = f(c);

        // Convergence check
        if (std::abs(fc) < criteria.abs_ftol) {
            return RootFindingResult<T>::success(c, iterations, std::abs(fc));
        }

        if (criteria.is_interval_converged(a, b)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(b - a));
        }

        // Bracket update + Anderson-Björk modification
        if (fa * fc < T(0)) {
            // The root is on the [a, c] side → b ← c
            b = c; fb = fc;
            if (side == 1) {
                // The a-side was not changed twice in a row → Anderson-Björk modification
                T m = T(1) - fc / fb;
                if (m <= T(0)) m = T(1) / T(2);
                fa *= m;
            }
            side = 1;
        } else {
            // The root is on the [c, b] side → a ← c
            a = c; fa = fc;
            if (side == -1) {
                T m = T(1) - fc / fa;
                if (m <= T(0)) m = T(1) / T(2);
                fb *= m;
            }
            side = -1;
        }
    }

    T best = (std::abs(fa) < std::abs(fb)) ? a : b;
    T best_f = std::min(std::abs(fa), std::abs(fb));
    return RootFindingResult<T>::partial_success(best, false, iterations, best_f);
}

// ---- Alefeld-Potra-Shi (TOMS 748) internal helpers ----
namespace detail {

/// Four-point inverse cubic interpolation (Alefeld-Potra-Shi)
/// Returns P(0) of the inverse cubic polynomial P(y) through (fa,a),(fb,b),(fd,d),(fe,e)
template<concepts::OrderedField T>
T aps_ipzero(T a, T b, T d, T e, T fa, T fb, T fd, T fe) {
    T q11 = (d - e) * fd / (fe - fd);
    T q21 = (b - d) * fb / (fd - fb);
    T q31 = (a - b) * fa / (fb - fa);
    T d21 = (b - d) * fd / (fd - fb);
    T d31 = (a - b) * fb / (fb - fa);
    T q22 = (d21 - q11) * fb / (fe - fb);
    T q32 = (d31 - q21) * fa / (fd - fa);
    T d32 = (d31 - q21) * fd / (fd - fa);
    T q33 = (d32 - q22) * fa / (fe - fa);
    return a + q31 + q32 + q33;
}

/// Find the root of the quadratic interpolation polynomial of f using k steps of Newton's method
template<concepts::OrderedField T>
T aps_newton_quadratic(T a, T b, T d, T fa, T fb, T fd, int k) {
    T B0 = (fb - fa) / (b - a);
    T C0 = ((fd - fb) / (d - b) - B0) / (d - a);
    T c = (std::abs(fa) <= std::abs(fb)) ? a : b;
    for (int i = 0; i < k; ++i) {
        T pc = fa + (c - a) * (B0 + (c - b) * C0);
        T pdc = B0 + C0 * (T(2) * c - a - b);
        if (std::abs(pdc) < std::numeric_limits<T>::epsilon()) break;
        c -= pc / pdc;
    }
    return c;
}

} // namespace detail

/**
 * @brief Root finding by the Alefeld-Potra-Shi method (TOMS 748)
 *
 * An improvement on Brent's method. Through inverse cubic interpolation + Newton
 * quadratic interpolation + bisection fallback, it achieves high-order convergence
 * while guaranteeing at least the convergence speed of bisection in the worst case.
 * Reference: Alefeld, Potra, Shi "Algorithm 748" (ACM TOMS, 1995)
 *
 * @param f The function whose root is sought
 * @param a Lower bound of the interval (requires f(a)·f(b) < 0)
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 */
template<concepts::OrderedField T>
RootFindingResult<T> alefeld_potra_shi(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a), fb = f(b);
    size_t nfe = 2;

    if (fa * fb > T(0))
        return RootFindingResult<T>::failure(0, std::abs(b - a));
    if (criteria.is_function_converged(fa))
        return RootFindingResult<T>::success(a, 0, std::numeric_limits<T>::epsilon());
    if (criteria.is_function_converged(fb))
        return RootFindingResult<T>::success(b, 0, std::numeric_limits<T>::epsilon());

    // Guarantee a < b
    if (a > b) { std::swap(a, b); std::swap(fa, fb); }

    T d{}, fd{}, e{}, fe{};
    T last_c{};

    // Restrict c to the interior of the interval (a, b) (keep it 2.5% inside the endpoints)
    auto clamp = [&](T c) -> T {
        T tol = (b - a) * T(0.025);
        c = std::max(a + tol, std::min(b - tol, c));
        if (c <= a || c >= b) c = (a + b) / T(2);
        return c;
    };

    // Whether all four function values are distinct (applicability condition for inverse cubic interpolation)
    auto use_cubic = [&]() -> bool {
        return fa != fb && fa != fd && fa != fe &&
               fb != fd && fb != fe && fd != fe;
    };

    // Return the best root candidate
    auto best = [&]() -> T {
        return (std::abs(fb) <= std::abs(fa)) ? b : a;
    };

    // Evaluate → check → bracket update. Return value: 0=continue, 1=f(c) converged, 2=interval converged
    auto step = [&](T c, bool update_e) -> int {
        c = clamp(c);
        last_c = c;
        T fc = f(c); ++nfe;
        if (criteria.is_function_converged(fc)) return 1;
        if (update_e) { e = d; fe = fd; }
        if (fa * fc < T(0)) { d = b; fd = fb; b = c; fb = fc; }
        else                 { d = a; fd = fa; a = c; fa = fc; }
        return criteria.is_interval_converged(a, b) ? 2 : 0;
    };

    // --- Step 1: secant method (e does not exist yet) ---
    int r = step(a - fa * (b - a) / (fb - fa), false);
    if (r == 1) return RootFindingResult<T>::success(last_c, nfe, std::abs(b - a));
    if (r == 2) return RootFindingResult<T>::success(best(), nfe, std::abs(b - a));

    // --- Step 2: Newton quadratic interpolation (set e for the first time) ---
    r = step(detail::aps_newton_quadratic(a, b, d, fa, fb, fd, 2), true);
    if (r == 1) return RootFindingResult<T>::success(last_c, nfe, std::abs(b - a));
    if (r == 2) return RootFindingResult<T>::success(best(), nfe, std::abs(b - a));

    // --- Main loop: four points available ---
    while (nfe < criteria.max_iterations) {
        T mu = (b - a) / T(2);

        // Three interpolation steps: (a) inverse cubic, (b) inverse cubic, (c) Newton quadratic
        for (int sub = 0; sub < 3; ++sub) {
            T c;
            if (sub < 2) {
                c = use_cubic()
                    ? detail::aps_ipzero(a, b, d, e, fa, fb, fd, fe)
                    : detail::aps_newton_quadratic(a, b, d, fa, fb, fd, 3);
            } else {
                c = detail::aps_newton_quadratic(a, b, d, fa, fb, fd, 2);
            }
            r = step(c, true);
            if (r == 1) return RootFindingResult<T>::success(last_c, nfe, std::abs(b - a));
            if (r == 2) return RootFindingResult<T>::success(best(), nfe, std::abs(b - a));
        }

        // If the interval did not halve, shrink it safely by bisection
        if ((b - a) >= mu) {
            r = step((a + b) / T(2), true);
            if (r == 1) return RootFindingResult<T>::success(last_c, nfe, std::abs(b - a));
            if (r == 2) return RootFindingResult<T>::success(best(), nfe, std::abs(b - a));
        }
    }

    return RootFindingResult<T>::partial_success(best(), false, nfe, std::abs(b - a));
}

// Wrapper functions for compatibility with the old interface (legacy support)
template<concepts::OrderedField T>
std::optional<T> bisection_legacy(
    const std::function<T(T)>& f, T a, T b,
    T tolerance = std::numeric_limits<T>::epsilon() * 100,
    size_t max_iterations = 100) {

    ConvergenceCriteria<T> criteria;
    criteria.abs_ftol = tolerance;
    criteria.abs_xtol = tolerance;
    criteria.max_iterations = max_iterations;

    auto result = bisection(f, a, b, criteria);
    return result.root;
}

template<concepts::OrderedField T>
std::optional<T> newton_raphson_legacy(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    T tolerance = std::numeric_limits<T>::epsilon() * 100,
    size_t max_iterations = 50) {

    ConvergenceCriteria<T> criteria;
    criteria.abs_ftol = tolerance;
    criteria.abs_xtol = tolerance;
    criteria.max_iterations = max_iterations;

    auto result = newton_raphson(f, df, x0, criteria);
    return result.root;
}

template<concepts::OrderedField T>
std::optional<T> secant_method_legacy(
    const std::function<T(T)>& f,
    T x0, T x1,
    T tolerance = std::numeric_limits<T>::epsilon() * 100,
    size_t max_iterations = 50) {

    ConvergenceCriteria<T> criteria;
    criteria.abs_ftol = tolerance;
    criteria.abs_xtol = tolerance;
    criteria.max_iterations = max_iterations;

    auto result = secant_method(f, x0, x1, criteria);
    return result.root;
}

template<concepts::OrderedField T>
std::optional<T> brent_method_legacy(
    const std::function<T(T)>& f,
    T a, T b,
    T tolerance = std::numeric_limits<T>::epsilon() * 100,
    size_t max_iterations = 100) {

    ConvergenceCriteria<T> criteria;
    criteria.abs_ftol = tolerance;
    criteria.abs_xtol = tolerance;
    criteria.max_iterations = max_iterations;

    auto result = brent_method(f, a, b, criteria);
    return result.root;
}

template<concepts::OrderedField T>
std::optional<T> halley_method_legacy(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    T tolerance = std::numeric_limits<T>::epsilon() * 100,
    size_t max_iterations = 50) {

    ConvergenceCriteria<T> criteria;
    criteria.abs_ftol = tolerance;
    criteria.abs_xtol = tolerance;
    criteria.max_iterations = max_iterations;

    auto result = halley_method(f, df, d2f, x0, criteria);
    return result.root;
}

template<concepts::OrderedField T>
std::optional<T> schroder_method_legacy(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    T tolerance = std::numeric_limits<T>::epsilon() * 100,
    size_t max_iterations = 50) {

    ConvergenceCriteria<T> criteria;
    criteria.abs_ftol = tolerance;
    criteria.abs_xtol = tolerance;
    criteria.max_iterations = max_iterations;

    auto result = schroder_method(f, df, d2f, x0, criteria);
    return result.root;
}

// ============================================================================
// King's method (improved Pegasus method)
// ============================================================================

/**
 * @brief Root finding by King's method (improved Pegasus method)
 *
 * A kind of bracketing method. After finding the interpolation point with Regula
 * Falsi, on same-side stagnation it corrects the f value by f1 ← f1·f2/(f2+f) and
 * interpolates again to accelerate convergence. Super-linear convergence faster
 * than the Illinois method.
 *
 * Lineage: Illinois method (fixed weight 1/2) → Pegasus method (Dowell & Jarratt,
 * 1972; weight f2/(f2+f)) → King's method (1973). King's method is an improved
 * version that removes one or two of the slow sub-steps of Pegasus to raise the
 * asymptotic convergence efficiency (convergence order ~1.839 using only first-order
 * divided differences); the name king_method is not a misnomer.
 *
 * References:
 *   R. F. King, "An improved Pegasus method for root finding", BIT 13, 423-427 (1973)
 *   M. Dowell & P. Jarratt, "The Pegasus method for computing the root of an equation", BIT 12, 503-508 (1972)
 *
 * @param f The function whose root is sought
 * @param a Lower bound of the interval (requires f(a)·f(b) < 0)
 * @param b Upper bound of the interval
 * @param criteria Convergence criteria
 * @return The approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> king_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T fa = f(a);
    T fb = f(b);
    size_t iterations = 0;

    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(iterations, std::abs(b - a));
    }

    if (std::abs(fa) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(a, iterations, std::numeric_limits<T>::epsilon());
    }
    if (std::abs(fb) < criteria.abs_ftol) {
        return RootFindingResult<T>::success(b, iterations, std::numeric_limits<T>::epsilon());
    }

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        // f(a) and f(b) are nearly equal → fall back to the midpoint
        if (std::abs(fb - fa) < criteria.abs_ftol) {
            a = (a + b) / 2;
            fa = f(a);
            continue;
        }

        // Approximate root by linear interpolation
        T x = (a * fb - b * fa) / (fb - fa);
        T fx = f(x);

        // Convergence check
        if (std::abs(fx) < criteria.abs_ftol || std::abs(b - a) < criteria.abs_xtol) {
            return RootFindingResult<T>::success(x, iterations, std::abs(fx));
        }

        if (fb * fx < 0) {
            // The root is on the [x, b] side — normal update
            a = x;
            fa = fx;
        } else {
            // The root is on the [a, x] side — King modification: shrink f(a) to prevent stagnation
            // Re-interpolation loop
            while (true) {
                fa = fa * fb / (fb + fx);
                b = x;
                fb = fx;

                if (std::abs(fb - fa) < criteria.abs_ftol) break;

                x = (a * fb - b * fa) / (fb - fa);
                fx = f(x);

                if (std::abs(fx) < criteria.abs_ftol || std::abs(b - a) < criteria.abs_xtol) {
                    return RootFindingResult<T>::success(x, iterations, std::abs(fx));
                }

                if (fb * fx < 0) {
                    a = b;
                    fa = fb;
                    b = x;
                    fb = fx;
                    break;
                }
            }
        }
    }

    T best = (std::abs(fa) < std::abs(fb)) ? a : b;
    T best_f = std::min(std::abs(fa), std::abs(fb));
    return RootFindingResult<T>::partial_success(best, false, iterations, best_f);
}

// ============================================================================
// Popovski's method (seventh-order convergence)
// ============================================================================

/**
 * @brief Root finding by Popovski's method (seventh-order convergence)
 *
 * A variation of Neta's method. It achieves seventh-order convergence per iteration
 * through three stages: Newton step → two-point interpolation → inverse interpolation.
 * Both f(x) and f'(x) are required.
 *
 * Reference: McNamee & Pan, "Numerical Methods for Roots of Polynomials" Part 2, p289
 *
 * @param f The function whose root is sought
 * @param df Derivative of f
 * @param x0 Initial approximation
 * @param criteria Convergence criteria
 * @return The approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> popovski_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    T fx = f(x);

    for (size_t iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        if (std::abs(fx) < criteria.abs_ftol) {
            return RootFindingResult<T>::success(x, iterations, std::abs(fx));
        }

        T dfx = df(x);
        if (std::abs(dfx) < std::numeric_limits<T>::epsilon()) {
            return RootFindingResult<T>::partial_success(x, false, iterations, std::abs(fx));
        }

        // Step 1: Newton step → w
        T w = x - fx / dfx;
        T fw = f(w);

        // Step 2: two-point interpolation → z
        T denom = T(2) * fw - fx;
        if (std::abs(denom) < std::numeric_limits<T>::epsilon()) {
            x = w;
            fx = fw;
            continue;
        }
        T z = w + fw * (x - w) / denom;
        T fz = f(z);

        // Step 3: inverse interpolation → x_{n+1}
        T zw = z - w;
        T xw = x - w;
        T d1 = (fz - fw) * xw * fx;
        T d2 = (fw - fx) * zw * fz;
        T denom2 = d1 + d2;

        if (std::abs(denom2) < std::numeric_limits<T>::epsilon()) {
            x = z;
            fx = fz;
            continue;
        }

        T x_new = z - (zw - xw) * (fw - fx) * zw * fz / denom2;
        T fx_new = f(x_new);

        // Convergence check
        if (std::abs(fx_new) < criteria.abs_ftol ||
            std::abs(x_new - x) < criteria.abs_xtol) {
            return RootFindingResult<T>::success(x_new, iterations, std::abs(fx_new));
        }

        x = x_new;
        fx = fx_new;
    }

    return RootFindingResult<T>::partial_success(x, false, criteria.max_iterations, std::abs(fx));
}

// ============================================================================
// Swift-Lindfield bracket search
// ============================================================================

/**
 * @brief Bracket search by the Swift-Lindfield method
 *
 * Starting from an initial point a and step size h, it searches for a point b such
 * that f(a)·f(b) ≤ 0. If it moves in the direction where |f| increases, it reverses
 * direction and doubles the step size. Used as preprocessing for bracketing methods
 * (bisection, Brent, etc.).
 *
 * Reference: McNamee & Pan, "Numerical Methods for Roots of Polynomials" Part 2, p5
 *
 * @param f The function whose root is sought
 * @param a Starting point of the search
 * @param h Initial step size (may be positive or negative)
 * @param max_iterations Maximum number of iterations
 * @return The bracketing interval {a, b} (f(a)·f(b) ≤ 0) or nullopt
 */
template<concepts::OrderedField T>
struct BracketResult {
    T a;    ///< One end of the interval
    T b;    ///< The other end of the interval (f(a)·f(b) ≤ 0)
    size_t iterations;
};

template<concepts::OrderedField T>
std::optional<BracketResult<T>> swift_lindfield_bracket(
    const std::function<T(T)>& f,
    T a, T h,
    size_t max_iterations = 50) {

    T fa = f(a);
    T b = a;
    T fb;

    for (size_t i = 0; i < max_iterations; ++i) {
        b += h;
        fb = f(b);

        if (fa * fb <= 0) {
            return BracketResult<T>{a, b, i + 1};
        }

        if (std::abs(fa) < std::abs(fb)) {
            // f(b) increased → reverse direction (step size ×2)
            h *= T(-2);
        } else {
            // f(b) decreased → same direction with step size ×2
            h *= T(2);
            a = b;
            fa = fb;
        }
    }

    return std::nullopt;
}

/**
 * @brief Root finding by the fixed-point iteration method (Picard iteration)
 *
 * For an iteration map g, it repeats x_{n+1} = g(x_n) to find a fixed point of g
 * (g(x*) = x*, i.e. a solution of x = g(x)). Note that the argument g is the
 * "iteration function itself", not the equation "f(x)=0 that gives the root".
 *
 * If g is a contraction map in a neighborhood of x* (|g'(x*)| < 1), it converges
 * from any initial value in that neighborhood. When |g'(x*)| > 1 it diverges,
 * reaching the maximum number of iterations and returning as non-converged
 * (converged == false).
 *
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param g The iteration map (searches for the fixed point x = g(x); not the f of f(x)=0)
 * @param x0 Initial value
 * @param criteria Convergence criteria (uses |x_{n+1} - x_n|)
 * @return Result object containing the approximate fixed point and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> fixed_point_iteration(
    const std::function<T(T)>& g,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    size_t iterations = 0;

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        T x_next = g(x);

        // Convergence check by |x_{n+1} - x_n|
        if (criteria.is_variable_converged(x_next, x)) {
            return RootFindingResult<T>::success(x_next, iterations, std::abs(x_next - x));
        }

        x = x_next;
    }

    // Maximum number of iterations reached (diverging or slow convergence). Return the best approximation as non-converged.
    return RootFindingResult<T>::partial_success(x, false, iterations, std::abs(g(x) - x));
}

/**
 * @brief Root finding by Steffensen's method (derivative-free, quadratic convergence)
 *
 * Achieves the same quadratic convergence as Newton's method without using the
 * derivative. It applies Aitken's Δ² acceleration to fixed-point iteration, and
 * unlike the secant method does not need to keep two points.
 * Update formula: x_{n+1} = x_n − f(x_n)² / (f(x_n + f(x_n)) − f(x_n)).
 * The denominator f(x+f) − f(x) approximates f(x)·f'(x). Each iteration needs two function evaluations.
 *
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param x0 Initial value (start near the root; it may diverge if far away)
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> steffensen_method(
    const std::function<T(T)>& f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    size_t iterations = 0;

    for (iterations = 0; iterations < criteria.max_iterations; ++iterations) {
        T fx = f(x);
        if (criteria.is_function_converged(fx)) {
            return RootFindingResult<T>::success(x, iterations, std::abs(fx));
        }

        // Denominator ≈ f(x)·f'(x)
        T denom = f(x + fx) - fx;
        if (std::abs(denom) < std::numeric_limits<T>::epsilon() * T(10)) {
            if (criteria.is_function_converged(fx)) {
                return RootFindingResult<T>::success(x, iterations, std::abs(fx));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx));
        }

        T x_new = x - fx * fx / denom;

        if (criteria.is_variable_converged(x_new, x)) {
            return RootFindingResult<T>::success(x_new, iterations + 1, std::abs(x_new - x));
        }
        x = x_new;
    }

    return RootFindingResult<T>::failure(iterations, std::abs(f(x)));
}

/**
 * @brief Root finding by inverse quadratic interpolation (standalone version)
 *
 * It evaluates at f = 0 the quadratic interpolation polynomial of the "inverse
 * function" x = g(f) through the three points (x0,f0),(x1,f1),(x2,f2) to obtain the
 * next approximation. This makes the internal element of Brent's method usable
 * independently. Derivative-free with super-linear convergence (≈1.84 order), but
 * it becomes unstable when function values get close, so its proper use is in
 * combination with a safety net (Brent's method).
 *
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param x0,x1,x2 Three distinct initial points (the f values must also be distinct)
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> inverse_quadratic_interpolation(
    const std::function<T(T)>& f,
    T x0, T x1, T x2,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T f0 = f(x0), f1 = f(x1), f2 = f(x2);
    size_t iterations = 0;

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        if (criteria.is_function_converged(f2)) {
            return RootFindingResult<T>::success(x2, iterations, std::abs(f2));
        }

        T d01 = f0 - f1, d02 = f0 - f2, d12 = f1 - f2;
        const T tiny = std::numeric_limits<T>::epsilon() * T(10);
        if (std::abs(d01) < tiny || std::abs(d02) < tiny || std::abs(d12) < tiny) {
            // f values are close so inverse quadratic interpolation is not possible (a case that needs a safety net)
            if (criteria.is_function_converged(f2)) {
                return RootFindingResult<T>::success(x2, iterations, std::abs(f2));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(f2));
        }

        // Evaluate x = Σ x_i Π_{j≠i} f_j/(f_i − f_j) at f=0
        T x3 = x0 * f1 * f2 / (d01 * d02)
             + x1 * f0 * f2 / (-d01 * d12)
             + x2 * f0 * f1 / (d02 * d12);

        if (criteria.is_variable_converged(x3, x2)) {
            return RootFindingResult<T>::success(x3, iterations, std::abs(x3 - x2));
        }

        x0 = x1; f0 = f1;
        x1 = x2; f1 = f2;
        x2 = x3; f2 = f(x3);
    }

    return RootFindingResult<T>::failure(iterations, std::abs(f2));
}

/**
 * @brief Sidi's generalized secant method
 *
 * Using the derivative p'(x_last) at the most recent point of the Newton
 * interpolation polynomial p through the latest K+1 points, it updates
 * x_{n+1} = x_n − f(x_n)/p'(x_n). K=1 is the secant method; increasing K raises the
 * convergence order (approaching quadratic as K→∞). Derivative-free.
 *
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param x0,x1 Two initial points
 * @param criteria Convergence criteria
 * @param memory Number of points K+1 used for interpolation (default 3 = quadratic interpolation; 2 matches the secant method)
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> sidi_method(
    const std::function<T(T)>& f,
    T x0, T x1,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>(),
    size_t memory = 3) {

    if (memory < 2) memory = 2;
    std::vector<T> X = {x0, x1};
    std::vector<T> F = {f(x0), f(x1)};
    size_t iterations = 0;

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        const size_t K = X.size();
        T f_last = F[K - 1];
        if (criteria.is_function_converged(f_last)) {
            return RootFindingResult<T>::success(X[K - 1], iterations, std::abs(f_last));
        }

        // Newton divided differences dd[j] = f[x_0,...,x_j]
        std::vector<T> dd = F;
        for (size_t j = 1; j < K; ++j)
            for (size_t i = K - 1; i >= j; --i) {
                dd[i] = (dd[i] - dd[i - 1]) / (X[i] - X[i - j]);
                if (i == j) break;   // Avoid size_t underflow
            }

        // p'(x_last):  Σ_{j≥1} dd[j] · (Π_{i<j}(t−x_i))'|_{t=x_last}
        T xl = X[K - 1];
        T deriv = T(0);
        for (size_t j = 1; j < K; ++j) {
            T s = T(0);
            for (size_t m = 0; m < j; ++m) {
                T prod = T(1);
                for (size_t i = 0; i < j; ++i)
                    if (i != m) prod *= (xl - X[i]);
                s += prod;
            }
            deriv += dd[j] * s;
        }

        if (std::abs(deriv) < std::numeric_limits<T>::epsilon() * T(10)) {
            return RootFindingResult<T>::failure(iterations, std::abs(f_last));
        }

        T x_new = xl - f_last / deriv;
        if (criteria.is_variable_converged(x_new, xl)) {
            return RootFindingResult<T>::success(x_new, iterations, std::abs(x_new - xl));
        }

        // Sliding window: add the new point and drop the old one to keep K+1 points
        X.push_back(x_new);
        F.push_back(f(x_new));
        if (X.size() > memory) { X.erase(X.begin()); F.erase(F.begin()); }
    }

    return RootFindingResult<T>::failure(iterations, std::abs(F.back()));
}

/**
 * @brief Householder's method (order 3 = the order after Newton/Halley)
 *
 * The order d=3 member of the Householder family. An extension of Newton's method
 * (d=1) and Halley's method (d=2), with fourth-order convergence. Requires the
 * derivatives f', f'', f'''.
 *   x_{n+1} = x_n − (6 f f'² − 3 f² f'') / (6 f'³ − 6 f f' f'' + f² f''')
 *
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought
 * @param df,d2f,d3f First to third derivatives
 * @param x0 Initial value
 * @param criteria Convergence criteria
 * @return Result object containing the approximate root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<T> householder_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    const std::function<T(T)>& d3f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    T x = x0;
    size_t iterations = 0;

    for (iterations = 0; iterations < criteria.max_iterations; ++iterations) {
        T fx = f(x);
        if (criteria.is_function_converged(fx)) {
            return RootFindingResult<T>::success(x, iterations, std::abs(fx));
        }

        T d1 = df(x), d2 = d2f(x), d3 = d3f(x);
        T num = T(6) * fx * d1 * d1 - T(3) * fx * fx * d2;
        T den = T(6) * d1 * d1 * d1 - T(6) * fx * d1 * d2 + fx * fx * d3;

        if (std::abs(den) < std::numeric_limits<T>::epsilon() * T(10)) {
            if (criteria.is_function_converged(fx)) {
                return RootFindingResult<T>::success(x, iterations, std::abs(fx));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx));
        }

        T x_new = x - num / den;
        if (criteria.is_variable_converged(x_new, x)) {
            return RootFindingResult<T>::success(x_new, iterations + 1, std::abs(x_new - x));
        }
        x = x_new;
    }

    return RootFindingResult<T>::failure(iterations, std::abs(f(x)));
}

/**
 * @brief Müller's method (parabolic interpolation, complex-root capable)
 *
 * Among the roots of the parabola (quadratic polynomial) through three points, it
 * proceeds toward the one closer to the most recent point. Since it proceeds in
 * complex numbers even when the discriminant is negative, it can reach complex roots
 * even from real initial points (useful for finding complex roots of polynomials).
 * Convergence order ≈1.84, derivative-free. f must be evaluable with complex arguments.
 *
 * @tparam T A field type with an order structure (e.g. real numbers)
 * @param f The function whose root is sought (Complex<T> → Complex<T>)
 * @param z0,z1,z2 Three distinct initial points
 * @param criteria Convergence criteria (uses abs_ftol for |f| and abs_xtol for |Δz|)
 * @return Result object containing the approximate complex root and convergence status
 */
template<concepts::OrderedField T>
RootFindingResult<Complex<T>> muller_method(
    const std::function<Complex<T>(Complex<T>)>& f,
    Complex<T> z0, Complex<T> z1, Complex<T> z2,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>()) {

    using C = Complex<T>;
    using sangi::abs;
    using sangi::sqrt;
    const T tiny = std::numeric_limits<T>::epsilon() * T(10);

    C f0 = f(z0), f1 = f(z1), f2 = f(z2);
    size_t iterations = 0;

    for (iterations = 1; iterations <= criteria.max_iterations; ++iterations) {
        if (abs(f2) < criteria.abs_ftol) {
            return RootFindingResult<C>::success(z2, iterations, abs(f2));
        }

        C h0 = z1 - z0, h1 = z2 - z1;
        if (abs(h0) < tiny || abs(h1) < tiny) {
            return RootFindingResult<C>::partial_success(z2, abs(f2) < criteria.abs_ftol,
                                                         iterations, abs(f2));
        }

        C d0 = (f1 - f0) / h0;
        C d1 = (f2 - f1) / h1;
        C a = (d1 - d0) / (h1 + h0);
        C b = a * h1 + d1;
        C c = f2;

        C disc = sqrt(b * b - C(T(4)) * a * c);
        C den1 = b + disc, den2 = b - disc;
        C den = (abs(den1) >= abs(den2)) ? den1 : den2;   // Larger denominator = closer root
        if (abs(den) < tiny) {
            return RootFindingResult<C>::partial_success(z2, false, iterations, abs(f2));
        }

        C z3 = z2 - (C(T(2)) * c) / den;

        if (abs(z3 - z2) < criteria.abs_xtol) {
            return RootFindingResult<C>::success(z3, iterations, abs(z3 - z2));
        }

        z0 = z1; f0 = f1;
        z1 = z2; f1 = f2;
        z2 = z3; f2 = f(z3);
    }

    return RootFindingResult<C>::partial_success(z2, false, iterations, abs(f2));
}

#if SANGI_HAS_MKL
// Add the MKL implementation here
#endif

} // namespace sangi

#endif // SANGI_ROOT_FINDING_1D_HPP