// 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"

namespace sangi {

/**
 * @brief 二分法による求根
 * @tparam T 順序構造を持つ体型（実数など）
 * @param f 根を求める関数
 * @param a 区間の下限
 * @param b 区間の上限
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態を含む結果オブジェクト
 */
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);

    // 区間[a, b]に根があるかチェック
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(0, std::abs(b - a));
    }

    // fa == 0 または fb == 0 の場合はその点が根
    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;  // 中点を計算
        T fc = f(c);

        // 複合的な収束判定
        if (criteria.is_function_converged(fc) || criteria.is_interval_converged(a, b)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(b - a) / 2);
        }

        // 区間を半分に絞り込む
        if (fa * fc < 0) {
            b = c;
            fb = fc;
        }
        else {
            a = c;
            fa = fc;
        }

        // 前回の中点からの変化がほとんどない場合も収束とみなす
        if (criteria.is_variable_converged(c, c_prev)) {
            return RootFindingResult<T>::success(c, iterations, std::abs(c - c_prev));
        }

        c_prev = c;
    }

    // 最大反復回数に達したが、最良の近似値を返す
    T c = (a + b) / 2;
    return RootFindingResult<T>::partial_success(c, false, iterations, std::abs(b - a) / 2);
}

/**
 * @brief ニュートン法による求根
 * @tparam T 順序構造を持つ体型（実数など）
 * @param f 根を求める関数
 * @param df 関数fの導関数
 * @param x0 初期値
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態を含む結果オブジェクト
 */
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);

        // 関数値による収束判定
        if (criteria.is_function_converged(fx, fx_prev)) {
            return RootFindingResult<T>::success(x, iterations, std::abs(x - x_prev));
        }

        T dfx = df(x);

        // 導関数がゼロに近い場合
        if (std::abs(dfx) < std::numeric_limits<T>::epsilon() * 10) {
            // 現在の関数値が十分小さければ収束とみなす
            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;

        // 変数値による収束判定
        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;
    }

    // 最大反復回数に達した場合
    return RootFindingResult<T>::failure(iterations, std::abs(f(x)));
}

/**
 * @brief 割線法による求根
 * @tparam T 順序構造を持つ体型（実数など）
 * @param f 根を求める関数
 * @param x0 初期値1
 * @param x1 初期値2
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態を含む結果オブジェクト
 */
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;

    // x0が既に根である場合
    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) {
        // 関数値による収束判定
        if (criteria.is_function_converged(fx1, fx0)) {
            return RootFindingResult<T>::success(x1, iterations, std::abs(x1 - x0));
        }

        // 割線の傾きがゼロに近い場合
        if (std::abs(fx1 - fx0) < std::numeric_limits<T>::epsilon()) {
            // 傾きがほぼゼロなら、fx1が十分小さければ収束とみなす
            if (criteria.is_function_converged(fx1)) {
                return RootFindingResult<T>::success(x1, iterations, std::abs(fx1));
            }
            return RootFindingResult<T>::failure(iterations, std::abs(fx1));
        }

        // 割線法の更新式
        T x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0);

        // 変数値による収束判定
        if (criteria.is_variable_converged(x2, x1)) {
            return RootFindingResult<T>::success(x2, iterations, std::abs(x2 - x1));
        }

        // 値の更新
        x0 = x1;
        fx0 = fx1;
        x1 = x2;
        fx1 = f(x1);
    }

    // 最大反復回数に達したが、収束していない場合
    return RootFindingResult<T>::failure(iterations, std::abs(fx1));
}

/**
 * @brief ブレント法による求根
 * @tparam T 順序構造を持つ体型（実数など）
 * @param f 根を求める関数
 * @param a 区間の下限
 * @param b 区間の上限
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態を含む結果オブジェクト
 */
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;

    // 区間[a, b]に根があるかチェック
    if (fa * fb > 0) {
        return RootFindingResult<T>::failure(iterations, std::abs(b - a));
    }

    // fa == 0 または 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());
    }

    // |fa| > |fb| となるように a と b を入れ替える
    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;  // 初期近似値はbから始める
    T fs = fb;
    T s_prev = a;  // 前回値は明確に異なる値に設定

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

        // 区間幅による収束判定
        if (criteria.is_interval_converged(a, b) && iterations > 1) {
            return RootFindingResult<T>::success(s, iterations, std::abs(b - a));
        }

        // 変数値による収束判定
        if (iterations > 1 && criteria.is_variable_converged(s, s_prev)) {
            return RootFindingResult<T>::success(s, iterations, std::abs(s - s_prev));
        }

        // 次の近似点を計算
        s_prev = s;
        T s_new;

        // 逆二次補間または割線法の計算
        if (std::abs(fc - fa) > criteria.abs_ftol &&
            std::abs(fb - fc) > criteria.abs_ftol) {
            // 逆二次補間法
            s_new = a * fb * fc / ((fa - fb) * (fa - fc)) +
                b * fa * fc / ((fb - fa) * (fb - fc)) +
                c * fa * fb / ((fc - fa) * (fc - fb));
        }
        else {
            // 割線法
            s_new = b - fb * (b - a) / (fb - fa);
        }

        // s_newが境界外、または良くない場合は二分法を使用
        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;  // 二分法
            mflag = true;
        }
        else {
            mflag = false;
        }

        // 新しい近似値での関数評価
        s = s_new;
        fs = f(s);

        // 点の更新
        d = c;
        c = b;
        fc = fb;

        // 根を含む新しい区間に更新
        if (fa * fs < 0) {
            b = s;
            fb = fs;
        }
        else {
            a = s;
            fa = fs;
        }

        // |fa| >= |fb| となるように点を入れ替え
        if (std::abs(fa) < std::abs(fb)) {
            std::swap(a, b);
            std::swap(fa, fb);
        }
    }

    // 最大反復回数に達した場合
    // 最も関数値が小さい点を返す
    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 Ridders 法による求根
 *
 * ブラケット法の一種。二分法の中点に指数関数変換を適用し、2次収束を実現する。
 * x_new = x_mid + (x_mid - a) * sign(fa - fb) * f_mid / sqrt(f_mid² - fa*fb)
 *
 * @param f 根を求める関数
 * @param a 区間の下限 (f(a)·f(b) < 0 が必要)
 * @param b 区間の上限
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態を含む結果オブジェクト
 */
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;

    // 区間[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 の公式
        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);

        // 収束判定
        if (std::abs(fx_new) < criteria.abs_ftol) {
            return RootFindingResult<T>::success(x_new, iterations, std::abs(fx_new));
        }

        // ブラケット更新
        if (fm * fx_new < T(0)) {
            // mid と 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;
        }

        // 区間幅による収束判定
        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 はさみうち法 (Illinois 改良型) による求根
 *
 * Regula Falsi (偽位法) の改良版。線形補間で次の近似点を求め、
 * 同側停滞時に Illinois 修正 (重み半減) を適用して収束を保証する。
 *
 * @param f 根を求める関数
 * @param a 区間の下限 (f(a)·f(b) < 0 が必要)
 * @param b 区間の上限
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態を含む結果オブジェクト
 */
template<concepts::OrderedField T>
RootFindingResult<T> false_position(
    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;

    // 区間[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=初期, 1=a側変化, -1=b側変化 (Illinois修正用)

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

        // 収束判定
        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));
        }

        // ブラケット更新 + Illinois 修正
        if (fa * fc < T(0)) {
            b = c; fb = fc;
            if (side == 1) {
                // a 側が2回連続で変更されなかった → Illinois 修正
                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 Halley 法による求根 (3次収束)
 *
 * f, f', f'' を用いて Newton 法より高次の収束を実現する。
 * 更新式: x_{n+1} = x_n - 2·f·f' / (2·f'² - f·f'')
 *
 * @param f   根を求める関数
 * @param df  f の1階導関数
 * @param d2f f の2階導関数
 * @param x0  初期値
 * @param 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);

        // 分母: 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 Schroder 法による求根 (3次収束, Householder 2次)
 *
 * f, f', f'' を用いる Householder 法の2次形態。
 * 重根に対しても2次収束を維持する。
 * 更新式: x_{n+1} = x_n - f·f' / (f'² - f·f'')
 *
 * @param f   根を求める関数
 * @param df  f の1階導関数
 * @param d2f f の2階導関数
 * @param x0  初期値
 * @param 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);

        // 分母: 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)));
}

// ---- TOMS 748 内部ヘルパー ----
namespace detail {

/// 4点逆3次補間 (Alefeld-Potra-Shi)
/// (fa,a),(fb,b),(fd,d),(fe,e) を通る逆3次多項式 P(y) の P(0) を返す
template<concepts::OrderedField T>
T toms748_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;
}

/// f の2次補間多項式の根を Newton 法 k ステップで求める
template<concepts::OrderedField T>
T toms748_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 TOMS 748 法 (Alefeld-Potra-Shi) による求根
 *
 * Brent 法の改良。逆3次補間 + Newton 2次補間 + 二分法フォールバックにより、
 * 高次収束を実現しつつ worst case でも二分法以上の収束速度を保証する。
 * 参考: Alefeld, Potra, Shi "Algorithm 748" (ACM TOMS, 1995)
 *
 * @param f 根を求める関数
 * @param a 区間の下限 (f(a)·f(b) < 0 が必要)
 * @param b 区間の上限
 * @param criteria 収束判定基準
 */
template<concepts::OrderedField T>
RootFindingResult<T> toms748(
    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());

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

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

    // c を区間 (a, b) の内部に制限 (端点の 2.5% 内側に留める)
    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;
    };

    // 4つの関数値がすべて異なるか (逆3次補間の適用条件)
    auto use_cubic = [&]() -> bool {
        return fa != fb && fa != fd && fa != fe &&
               fb != fd && fb != fe && fd != fe;
    };

    // 最良の根候補を返す
    auto best = [&]() -> T {
        return (std::abs(fb) <= std::abs(fa)) ? b : a;
    };

    // 評価 → 判定 → ブラケット更新。戻り値: 0=続行, 1=f(c)収束, 2=区間収束
    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: セカント法 (e はまだ存在しない) ---
    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 2次補間 (e を初めて設定) ---
    r = step(detail::toms748_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));

    // --- メインループ: 4点利用可能 ---
    while (nfe < criteria.max_iterations) {
        T mu = (b - a) / T(2);

        // 3回の補間ステップ: (a) 逆3次, (b) 逆3次, (c) Newton 2次
        for (int sub = 0; sub < 3; ++sub) {
            T c;
            if (sub < 2) {
                c = use_cubic()
                    ? detail::toms748_ipzero(a, b, d, e, fa, fb, fd, fe)
                    : detail::toms748_newton_quadratic(a, b, d, fa, fb, fd, 3);
            } else {
                c = detail::toms748_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 ((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));
}

// 旧インターフェースとの互換性のためのラッパー関数（レガシーサポート）
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 法 (改良 Pegasus 法)
// ============================================================================

/**
 * @brief King 法による求根 (改良 Pegasus 法)
 *
 * ブラケット法の一種。Regula Falsi で補間点を求めた後、同側停滞時に
 * f 値を f1 ← f1·f2/(f2+f) で修正し、再度補間を行うことで収束を加速する。
 * Illinois 法より速い超線形収束。
 *
 * 参考: R.F.King, "An improved Pegasus method for root finding", BIT 1973
 *
 * @param f 根を求める関数
 * @param a 区間の下限 (f(a)·f(b) < 0 が必要)
 * @param b 区間の上限
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態
 */
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) と f(b) がほぼ等しい → 中点に退避
        if (std::abs(fb - fa) < criteria.abs_ftol) {
            a = (a + b) / 2;
            fa = f(a);
            continue;
        }

        // 線形補間で近似根
        T x = (a * fb - b * fa) / (fb - fa);
        T 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) {
            // 根は [x, b] 側 — 通常の更新
            a = x;
            fa = fx;
        } else {
            // 根は [a, x] 側 — King 修正: f(a) を縮小して停滞を防ぐ
            // 再補間ループ
            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 法 (7次収束)
// ============================================================================

/**
 * @brief Popovski 法による求根 (7次収束)
 *
 * Neta 法のバリエーション。Newton ステップ → 2点補間 → 逆補間の3段階で
 * 各反復 7次収束を達成。f(x) と f'(x) の両方が必要。
 *
 * 参考: McNamee & Pan, "Numerical Methods for Roots of Polynomials" Part 2, p289
 *
 * @param f 根を求める関数
 * @param df f の導関数
 * @param x0 初期近似値
 * @param criteria 収束判定基準
 * @return 根の近似値と収束状態
 */
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 ステップ → w
        T w = x - fx / dfx;
        T fw = f(w);

        // Step 2: 2点補間 → 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: 逆補間 → 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);

        // 収束判定
        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 ブラケット探索
// ============================================================================

/**
 * @brief Swift-Lindfield 法によるブラケット探索
 *
 * 初期点 a と刻み幅 h から出発し、f(a)·f(b) ≤ 0 となる点 b を探索する。
 * |f| が増大する方向に進んだ場合は逆方向に転換し刻み幅を倍増する。
 * ブラケット法 (bisection, Brent 等) の前処理として使用。
 *
 * 参考: McNamee & Pan, "Numerical Methods for Roots of Polynomials" Part 2, p5
 *
 * @param f 根を求める関数
 * @param a 探索の開始点
 * @param h 初期刻み幅 (正でも負でもよい)
 * @param max_iterations 最大反復回数
 * @return ブラケット区間 {a, b} (f(a)·f(b) ≤ 0) または nullopt
 */
template<concepts::OrderedField T>
struct BracketResult {
    T a;    ///< 区間の一端
    T b;    ///< 区間の他端 (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) が大きくなった → 逆方向に転換 (刻み幅 2倍)
            h *= T(-2);
        } else {
            // f(b) が小さくなった → 同方向に刻み幅 2倍
            h *= T(2);
            a = b;
            fa = fb;
        }
    }

    return std::nullopt;
}

#if SANGI_HAS_MKL
// MKL実装をここに追加
#endif

} // namespace sangi

#endif // SANGI_ROOT_FINDING_1D_HPP