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

// polynomial_roots.hpp
// 多項式の数値的求根アルゴリズム
// lib++20 代数方程式/* からの移植 (2026-02-18)
//
// 低次 (1〜4次): 桁落ちに対して頑健な公式解
// 高次 (5次以上): Jenkins-Traub, Laguerre, DKA 反復法
//
// 全関数は Complex<T> のベクトルを返す (実根は im=0)

#ifndef SANGI_POLYNOMIAL_ROOTS_HPP
#define SANGI_POLYNOMIAL_ROOTS_HPP

#include <math/core/Polynomial.hpp>
#include <math/core/Complex.hpp>

#include <vector>
#include <cmath>
#include <algorithm>
#include <limits>
#include <cassert>

namespace sangi {

// ================================================================
// 内部ヘルパー
// ================================================================

namespace detail {

// 符号関数
template<typename T>
T sgn(T val) {
    if (val > T(0)) return T(1);
    if (val < T(0)) return T(-1);
    return T(0);
}

// 立方根 (実数)
template<typename T>
T cubeRoot(T x) {
    using std::cbrt;
    return cbrt(x);
}

} // namespace detail

// ================================================================
// A-1. solveLinear (1次方程式)
// ================================================================

/// 1次方程式 a*x + b = 0 の解を返す
/// p.coefficients(): [b, a] (昇べき)
template<typename T>
std::vector<Complex<T>> solveLinear(const Polynomial<T>& p) {
    assert(p.degree() == 1);

    T a = p[1];  // x の係数
    T b = p[0];  // 定数項

    T root = -b / a;
    return { Complex<T>(root) };
}

// ================================================================
// A-2. solveQuadratic (2次方程式)
// ================================================================

/// 2次方程式 a*x^2 + b*x + c = 0 の解を返す
/// 桁落ち対策: 大きい絶対値の根を先に求め、Vieta の公式で小さい根を計算
/// 【出典】奥村「アルゴリズム事典」p205
template<typename T>
std::vector<Complex<T>> solveQuadratic(const Polynomial<T>& p) {
    assert(p.degree() == 2);

    T a = p[2];  // x^2 の係数
    T b = p[1];  // x の係数
    T c = p[0];  // 定数項

    // モニック化: x^2 + px + q = 0
    T pp = b / a;
    T q  = c / a;

    using std::sqrt;
    using std::abs;

    if (q == T(0)) {
        // x^2 + px = x(x+p) = 0
        return { Complex<T>(T(0)), Complex<T>(-pp) };
    }

    T halfP = pp / T(2);          // p/2
    T disc  = halfP * halfP - q;  // 判別式 D/4

    if (disc > T(0)) {
        // 実根2個 — 桁落ち対策
        // 絶対値が大きい方の根を先に計算
        T x1;
        if (halfP > T(0))
            x1 = -halfP - sqrt(disc);
        else
            x1 = -halfP + sqrt(disc);

        T x2 = q / x1;  // Vieta: x1*x2 = q
        return { Complex<T>(x1), Complex<T>(x2) };
    }
    else if (disc < T(0)) {
        // 複素共役根
        T re = -halfP;
        T im = sqrt(-disc);
        return { Complex<T>(re, im), Complex<T>(re, -im) };
    }
    else {
        // 重根
        T x = -halfP;
        return { Complex<T>(x), Complex<T>(x) };
    }
}

// ================================================================
// A-3. solveCubic (3次方程式)
// ================================================================

/// 3次方程式 a*x^3 + b*x^2 + c*x + d = 0 の解を返す
/// Cardano の公式 + 三角関数法 (casus irreducibilis)
template<typename T>
std::vector<Complex<T>> solveCubic(const Polynomial<T>& p) {
    assert(p.degree() == 3);

    using std::abs;
    using std::sqrt;
    using std::cbrt;
    using std::cos;
    using std::acos;
    using std::pow;

    const T pi = std::acos(T(-1));

    T a = p[3];
    T b = p[2];
    T c = p[1];
    T d = p[0];

    // モニック化
    b /= a;  c /= a;  d /= a;

    // Tschirnhaus 変換: t = x + b/3 → t^3 + pt + q = 0
    T b3 = b / T(3);
    T pp = c - b * b3;                           // p = c - b^2/3
    T qq = d - b3 * c + T(2) * b3 * b3 * b3;    // q = d - bc/3 + 2b^3/27

    // 判別式 Δ = -(4p^3 + 27q^2)
    T disc = -(T(4) * pp * pp * pp + T(27) * qq * qq);

    std::vector<Complex<T>> roots;

    if (disc > T(0)) {
        // Casus irreducibilis: 実根3個 → 三角関数法
        T r = sqrt(-pp / T(3));                        // r = √(-p/3)
        T cosTheta = -qq / (T(2) * r * r * r);        // cos(θ) = -q/(2r^3)
        // 数値誤差で ±1 を超える場合をクランプ
        if (cosTheta > T(1))  cosTheta = T(1);
        if (cosTheta < T(-1)) cosTheta = T(-1);
        T theta = acos(cosTheta);

        T t0 = T(2) * r * cos(theta / T(3));
        T t1 = T(2) * r * cos((theta + T(2) * pi) / T(3));
        T t2 = T(2) * r * cos((theta + T(4) * pi) / T(3));

        roots.push_back(Complex<T>(t0 - b3));
        roots.push_back(Complex<T>(t1 - b3));
        roots.push_back(Complex<T>(t2 - b3));
    }
    else if (disc < T(0)) {
        // 実根1個 + 複素共役根
        T sqrtQ = qq * qq / T(4) + pp * pp * pp / T(27);
        T sqrtVal = sqrt(sqrtQ);

        T A = -qq / T(2) + sqrtVal;
        T B = -qq / T(2) - sqrtVal;
        A = detail::cubeRoot(A);
        B = detail::cubeRoot(B);

        T realRoot = A + B - b3;
        roots.push_back(Complex<T>(realRoot));

        // 残り2根: (x - realRoot) で割って2次式を解く
        // t^3 + pt + q = (t - (A+B))(t^2 + (A+B)t + ...) から
        // 2次式の係数を合成除算で求める
        T sumAB = A + B;
        T re = -sumAB / T(2) - b3;
        T im = abs(A - B) * sqrt(T(3)) / T(2);

        roots.push_back(Complex<T>(re, im));
        roots.push_back(Complex<T>(re, -im));
    }
    else {
        // Δ = 0: 重根あり
        if (pp == T(0) && qq == T(0)) {
            // 三重根: t^3 = 0
            T x = -b3;
            roots.push_back(Complex<T>(x));
            roots.push_back(Complex<T>(x));
            roots.push_back(Complex<T>(x));
        }
        else {
            // 二重根 + 単根
            // t³+pt+q=0 が (t-α)²(t+2α)=0 のとき p=-3α², q=2α³
            // → α = -3q/(2p), β = -2α = 3q/p
            T doubleRoot = -T(3) * qq / (T(2) * pp);
            T singleRoot = T(3) * qq / pp;
            roots.push_back(Complex<T>(singleRoot - b3));
            roots.push_back(Complex<T>(doubleRoot - b3));
            roots.push_back(Complex<T>(doubleRoot - b3));
        }
    }

    return roots;
}

// ================================================================
// A-4. solveQuartic (4次方程式)
// ================================================================

/// 4次方程式 a*x^4 + b*x^3 + c*x^2 + d*x + e = 0 の解を返す
/// Ferrari の方法: 解決3次式の実根を用いて2つの2次式に分解
template<typename T>
std::vector<Complex<T>> solveQuartic(const Polynomial<T>& p) {
    assert(p.degree() == 4);

    using std::abs;
    using std::sqrt;

    T a4 = p[4];
    T a3 = p[3];
    T a2 = p[2];
    T a1 = p[1];
    T a0 = p[0];

    // モニック化
    a3 /= a4;  a2 /= a4;  a1 /= a4;  a0 /= a4;

    // Depressed quartic: t = x - a3/4
    // t^4 + pt^2 + qt + r = 0
    T shift = a3 / T(4);
    T pp = a2 - T(6) * shift * shift;
    T qq = a1 - T(2) * a2 * shift + T(8) * shift * shift * shift;
    T rr = a0 - a1 * shift + a2 * shift * shift - T(3) * shift * shift * shift * shift;

    std::vector<Complex<T>> roots;

    if (abs(qq) < std::numeric_limits<T>::epsilon() * T(100)) {
        // qq ≈ 0 → biquadratic: t^4 + pt^2 + r = 0
        // u = t^2 として u^2 + pu + r = 0
        Polynomial<T> biQuad(std::vector<T>{rr, pp, T(1)});
        auto uRoots = solveQuadratic(biQuad);

        for (const auto& u : uRoots) {
            auto sqrtU = sangi::sqrt(u);
            roots.push_back(Complex<T>(sqrtU.re - shift, sqrtU.im));
            roots.push_back(Complex<T>(-sqrtU.re - shift, -sqrtU.im));
        }
        return roots;
    }

    // 解決3次式 (resolvent cubic): 8y^3 - 4py^2 - 8ry + (4pr - q^2) = 0
    Polynomial<T> resolvent(std::vector<T>{
        T(4) * pp * rr - qq * qq,
        T(-8) * rr,
        T(-4) * pp,
        T(8)
    });

    auto cubicRoots = solveCubic(resolvent);

    // 実根の中で最大のものを選ぶ (数値安定性)
    T y0 = cubicRoots[0].re;
    for (size_t i = 1; i < cubicRoots.size(); ++i) {
        if (abs(cubicRoots[i].im) < std::numeric_limits<T>::epsilon() * T(100)) {
            if (cubicRoots[i].re > y0)
                y0 = cubicRoots[i].re;
        }
    }

    // 2つの2次式に分解
    // t^4 + pt^2 + qt + r = (t^2 + αt + β)(t^2 - αt + γ)
    // α^2 = 2y0 - p, β = y0 - q/(2α), γ = y0 + q/(2α)
    T alpha2 = T(2) * y0 - pp;
    T alpha;
    if (alpha2 > T(0)) {
        alpha = sqrt(alpha2);
    }
    else if (alpha2 > -std::numeric_limits<T>::epsilon() * T(100)) {
        alpha = T(0);
    }
    else {
        // フォールバック: 数値誤差による微小な負の値
        alpha = sqrt(abs(alpha2));
    }

    T beta, gamma;
    if (abs(alpha) > std::numeric_limits<T>::epsilon() * T(100)) {
        beta  = y0 - qq / (T(2) * alpha);
        gamma = y0 + qq / (T(2) * alpha);
    }
    else {
        // α ≈ 0 の場合
        T sqrtY0sq_r = sqrt(abs(y0 * y0 - rr));
        beta  = y0 - sqrtY0sq_r;
        gamma = y0 + sqrtY0sq_r;
    }

    // 2つの2次方程式を解く
    Polynomial<T> quad1(std::vector<T>{beta, alpha, T(1)});
    Polynomial<T> quad2(std::vector<T>{gamma, -alpha, T(1)});

    auto roots1 = solveQuadratic(quad1);
    auto roots2 = solveQuadratic(quad2);

    // シフトを戻す
    for (auto& r : roots1)
        roots.push_back(Complex<T>(r.re - shift, r.im));
    for (auto& r : roots2)
        roots.push_back(Complex<T>(r.re - shift, r.im));

    return roots;
}

// ================================================================
// B-1. Jenkins-Traub 法
// ================================================================

/// Jenkins-Traub 法による多項式求根 (TOMS493 ベース)
/// 最も頑健な多項式求根法。3段階アルゴリズム。
/// 移植元: lib++20 JenkinsTraub.cpp (C. Bond の実装ベース)
template<typename T>
std::vector<Complex<T>> jenkinsTraub(const Polynomial<T>& poly) {
    using std::abs;
    using std::sqrt;
    using std::log;
    using std::exp;
    using std::cos;
    using std::sin;

    int degree = poly.degree();
    if (degree < 1)
        return {};

    // 降べき順の係数配列を作成 (TOMS493 の慣習)
    // poly は昇べき: coeffs_[0]=定数, coeffs_[n]=主係数
    // op[0]=主係数, op[n]=定数項
    std::vector<T> op(degree + 1);
    for (int i = 0; i <= degree; ++i)
        op[i] = poly[degree - i];

    // 主係数がゼロなら不正
    if (op[0] == T(0))
        return {};

    // 結果格納用
    std::vector<T> zeror(degree), zeroi(degree);
    int foundCount = 0;

    // 作業用配列
    int n = degree;
    std::vector<T> p(n + 1), qp(n + 1), k(n + 1), qk(n + 1), svk(n + 1);
    std::vector<T> temp(n + 1), pt(n + 1);

    // 機械定数
    const T base = T(2);
    const T eta = std::numeric_limits<T>::epsilon();
    const T infin = std::numeric_limits<T>::max() / T(10);
    const T smalno = std::numeric_limits<T>::min() * T(10);
    const T are = eta;
    const T mre = eta;
    const T lo = smalno / eta;

    // 回転角 (94度)
    T rot = T(94) * T(0.017453293);
    T cosr = cos(rot);
    T sinr = sin(rot);
    T xx = sqrt(T(0.5));
    T yy = -xx;

    // 零根を先に除去
    while (n > 0 && op[n] == T(0)) {
        int j = degree - n;
        zeror[j] = T(0);
        zeroi[j] = T(0);
        foundCount++;
        n--;
    }
    if (n < 1) goto done;

    // 係数を p にコピー
    for (int i = 0; i <= n; ++i)
        p[i] = op[i];

    // メインの求根ループ
    {
        // ローカル変数
        T sr, si, u, v, a, b, c, d;
        T a1, a3, a7, ee, f, g, h;
        T szr, szi, lzr, lzi;

        // --- 内部関数をラムダで定義 ---

        // quadsd: p を (1, u, v) の2次式で合成除算
        auto quadsd = [&](int nn, T uu, T vv, const std::vector<T>& pp,
                          std::vector<T>& qq, T& aa, T& bb) {
            bb = pp[0];
            qq[0] = bb;
            aa = pp[1] - bb * uu;
            qq[1] = aa;
            for (int i = 2; i <= nn; ++i) {
                T cc = pp[i] - aa * uu - bb * vv;
                qq[i] = cc;
                bb = aa;
                aa = cc;
            }
        };

        // quad: 2次方程式 a*z^2 + b1*z + c = 0 の解
        auto quad = [&](T qa, T b1, T qc,
                        T& sr_out, T& si_out, T& lr_out, T& li_out) {
            if (qa == T(0)) {
                sr_out = (b1 != T(0)) ? -qc / b1 : T(0);
                lr_out = T(0);
                si_out = T(0);
                li_out = T(0);
                return;
            }
            if (qc == T(0)) {
                sr_out = T(0);
                lr_out = -b1 / qa;
                si_out = T(0);
                li_out = T(0);
                return;
            }
            // 判別式のオーバーフロー回避
            T bb = b1 / T(2);
            T dd, eee;
            if (abs(bb) < abs(qc)) {
                eee = (qc < T(0)) ? -qa : qa;
                eee = bb * (bb / abs(qc)) - eee;
                dd = sqrt(abs(eee)) * sqrt(abs(qc));
            } else {
                eee = T(1) - (qa / bb) * (qc / bb);
                dd = sqrt(abs(eee)) * abs(bb);
            }
            if (eee < T(0)) {
                // 複素共役根
                sr_out = -bb / qa;
                lr_out = sr_out;
                si_out = abs(dd / qa);
                li_out = -si_out;
            } else {
                // 実根
                if (bb >= T(0)) dd = -dd;
                lr_out = (-bb + dd) / qa;
                sr_out = (lr_out != T(0)) ? (qc / lr_out) / qa : T(0);
                si_out = T(0);
                li_out = T(0);
            }
        };

        // calcsc: スカラー量の計算
        int calcType;
        auto calcsc = [&]() {
            quadsd(n - 1, u, v, k, qk, c, d);
            if (abs(c) > abs(k[n - 1]) * T(100) * eta ||
                abs(d) > abs(k[n - 2]) * T(100) * eta) {
                // type 1 or 2
                if (abs(d) < abs(c)) {
                    calcType = 1;
                    ee = a / c;
                    f = d / c;
                    g = u * ee;
                    h = v * b;
                    a3 = a * ee + (h / c + g) * b;
                    a1 = b - a * (d / c);
                    a7 = a + g * d + h * f;
                } else {
                    calcType = 2;
                    ee = a / d;
                    f = c / d;
                    g = u * b;
                    h = v * b;
                    a3 = (a + g) * ee + h * (b / d);
                    a1 = b * f - a;
                    a7 = (f + u) * a + h;
                }
            } else {
                calcType = 3;
            }
        };

        // nextk: 次の K 多項式を計算
        auto nextk = [&]() {
            if (calcType == 3) {
                k[0] = T(0);
                k[1] = T(0);
                for (int i = 2; i < n; ++i)
                    k[i] = qk[i - 2];
                return;
            }
            T tmp = (calcType == 1) ? b : a;
            if (abs(a1) <= abs(tmp) * eta * T(10)) {
                k[0] = T(0);
                k[1] = -a7 * qp[0];
                for (int i = 2; i < n; ++i)
                    k[i] = a3 * qk[i - 2] - a7 * qp[i - 1];
                return;
            }
            a7 /= a1;
            a3 /= a1;
            k[0] = qp[0];
            k[1] = qp[1] - a7 * qp[0];
            for (int i = 2; i < n; ++i)
                k[i] = a3 * qk[i - 2] - a7 * qp[i - 1] + qp[i];
        };

        // newest: 新しい u, v の推定
        auto newest = [&](T& uu, T& vv) {
            if (calcType == 3) {
                uu = T(0);
                vv = T(0);
                return;
            }
            T a4_, a5_;
            if (calcType == 2) {
                a4_ = (a + g) * f + h;
                a5_ = (f + u) * c + v * d;
            } else {
                a4_ = a + u * b + h * f;
                a5_ = c + (u + v * f) * d;
            }
            T b1_ = -k[n - 1] / p[n];
            T b2_ = -(k[n - 2] + b1_ * p[n - 1]) / p[n];
            T c1_ = v * b2_ * a1;
            T c2_ = b1_ * a7;
            T c3_ = b1_ * b1_ * a3;
            T c4_ = c1_ - c2_ - c3_;
            T tmp_ = a5_ + b1_ * a4_ - c4_;
            if (tmp_ == T(0)) {
                uu = T(0);
                vv = T(0);
                return;
            }
            uu = u - (u * (c3_ + c2_) + v * (b1_ * a1 + b2_ * a7)) / tmp_;
            vv = v * (T(1) + c4_ / tmp_);
        };

        // realit: 実零点の反復
        auto realit = [&](T sss, int& nz, int& iflag) {
            nz = 0;
            T s = sss;
            iflag = 0;
            int j = 0;
            T omp = T(0), t = T(0);

            while (true) {
                T pv = p[0];
                qp[0] = pv;
                for (int i = 1; i <= n; ++i) {
                    pv = pv * s + p[i];
                    qp[i] = pv;
                }
                T mp = abs(pv);
                T ms = abs(s);
                T eee = (mre / (are + mre)) * abs(qp[0]);
                for (int i = 1; i <= n; ++i)
                    eee = eee * ms + abs(qp[i]);

                if (mp <= T(20) * ((are + mre) * eee - mre * mp)) {
                    nz = 1;
                    szr = s;
                    szi = T(0);
                    return;
                }
                j++;
                if (j > 10) return;
                if (j >= 2) {
                    if (abs(t) <= T(0.001) * abs(s - t) && mp >= omp) {
                        iflag = 1;
                        sss = s;
                        return;
                    }
                }
                omp = mp;
                // K 多項式の更新
                T kv = k[0];
                qk[0] = kv;
                for (int i = 1; i < n; ++i) {
                    kv = kv * s + k[i];
                    qk[i] = kv;
                }
                if (abs(kv) <= abs(k[n - 1]) * T(10) * eta) {
                    k[0] = T(0);
                    for (int i = 1; i < n; ++i)
                        k[i] = qk[i - 1];
                } else {
                    t = -pv / kv;
                    k[0] = qp[0];
                    for (int i = 1; i < n; ++i)
                        k[i] = t * qk[i - 1] + qp[i];
                }
                kv = k[0];
                for (int i = 1; i < n; ++i)
                    kv = kv * s + k[i];
                t = T(0);
                if (abs(kv) > abs(k[n - 1]) * T(10) * eta)
                    t = -pv / kv;
                s += t;
            }
        };

        // quadit: 2次因子の反復
        auto quadit = [&](T uu, T vv, int& nz) {
            nz = 0;
            int tried = 0;
            u = uu;
            v = vv;
            int j = 0;
            T omp = T(0), relstp = T(0);

            while (true) {
                quad(T(1), u, v, szr, szi, lzr, lzi);
                if (abs(abs(szr) - abs(lzr)) > T(0.01) * abs(lzr))
                    return;
                quadsd(n, u, v, p, qp, a, b);
                T mp = abs(a - szr * b) + abs(szi * b);
                // 丸め誤差の上界
                T zm = sqrt(abs(v));
                T eee = T(2) * abs(qp[0]);
                T t = -szr * b;
                for (int i = 1; i < n; ++i)
                    eee = eee * zm + abs(qp[i]);
                eee = eee * zm + abs(a + t);
                eee *= (T(5) * mre + T(4) * are);
                eee -= (T(5) * mre + T(2) * are) * (abs(a) + t) + abs(b) * zm; // HVE fix
                eee += T(2) * are * abs(t);

                if (mp <= T(20) * eee) {
                    nz = 2;
                    return;
                }
                j++;
                if (j > 20) return;
                if (j >= 2 && !(relstp > T(0.01) || mp < omp || tried)) {
                    if (relstp < eta) relstp = eta;
                    relstp = sqrt(relstp);
                    u = u - u * relstp;
                    v = v + v * relstp;
                    quadsd(n, u, v, p, qp, a, b);
                    for (int i = 0; i < 5; ++i) {
                        calcsc();
                        nextk();
                    }
                    tried = 1;
                    j = 0;
                }
                omp = mp;
                calcsc();
                nextk();
                calcsc();
                T ui, vi;
                newest(ui, vi);
                if (vi == T(0)) return;
                relstp = abs((vi - v) / vi);
                u = ui;
                v = vi;
            }
        };

        // fxshfr: 固定シフト反復
        auto fxshfr = [&](int l2, int& nz) {
            nz = 0;
            T betav = T(0.25), betas = T(0.25);
            T oss = sr, ovv = v;
            T otv = T(1), ots = T(1);
            quadsd(n, u, v, p, qp, a, b);
            calcsc();

            for (int j = 0; j < l2; ++j) {
                nextk();
                calcsc();
                T ui, vi;
                newest(ui, vi);
                T vv = vi;
                T ss = T(0);
                if (k[n - 1] != T(0)) ss = -p[n] / k[n - 1];
                T tv = T(1), ts = T(1);
                if (j > 0 && calcType != 3) {
                    if (vv != T(0)) tv = abs((vv - ovv) / vv);
                    if (ss != T(0)) ts = abs((ss - oss) / ss);
                    T tvv = (tv < otv) ? tv * otv : T(1);
                    T tss = (ts < ots) ? ts * ots : T(1);
                    int vpass = (tvv < betav) ? 1 : 0;
                    int spass = (tss < betas) ? 1 : 0;

                    if (spass || vpass) {
                        T svu = u, svv_ = v;
                        for (int i = 0; i < n; ++i)
                            svk[i] = k[i];
                        T s = ss;
                        int vtry = 0, stry = 0;

                        // 収束の速い方を試す
                        bool tryQuad = !(spass && (!vpass || tss < tvv));

                        if (tryQuad) goto try_quad;
                        goto try_real;

                    try_real:
                        {
                            int iflag = 0;
                            realit(s, nz, iflag);
                            if (nz > 0) return;
                            stry = 1;
                            betas *= T(0.25);
                            if (iflag != 0) {
                                ui = -(s + s);
                                vi = s * s;
                                goto try_quad;
                            }
                            goto restore;
                        }

                    try_quad:
                        quadit(ui, vi, nz);
                        if (nz > 0) return;
                        vtry = 1;
                        betav *= T(0.25);
                        if (!stry && spass) {
                            for (int i = 0; i < n; ++i)
                                k[i] = svk[i];
                            goto try_real;
                        }

                    restore:
                        u = svu;
                        v = svv_;
                        for (int i = 0; i < n; ++i)
                            k[i] = svk[i];
                        if (vpass && !vtry) goto try_quad;

                        quadsd(n, u, v, p, qp, a, b);
                        calcsc();
                    }
                }
                ovv = vv;
                oss = ss;
                otv = tv;
                ots = ts;
            }
        };

        // --- メインループ ---
        auto findRoots = [&]() {
            while (n > 2) {
                // 係数スケーリング
                T maxCoef = T(0), minCoef = infin;
                for (int i = 0; i <= n; ++i) {
                    T x = abs(p[i]);
                    if (x > maxCoef) maxCoef = x;
                    if (x != T(0) && x < minCoef) minCoef = x;
                }
                T sc = lo / minCoef;
                bool doScale = true;
                if (sc > T(1) && infin / sc < maxCoef) doScale = false;
                if (sc <= T(1)) {
                    if (maxCoef < T(10)) doScale = false;
                    if (sc == T(0)) sc = smalno;
                }
                if (doScale) {
                    int l = static_cast<int>(log(sc) / log(base) + T(0.5));
                    T factor = (l >= 0) ? std::pow(base, static_cast<T>(l))
                                        : std::pow(T(1) / base, static_cast<T>(-l));
                    if (factor != T(1)) {
                        for (int i = 0; i <= n; ++i)
                            p[i] *= factor;
                    }
                }

                // 根の下界を計算
                for (int i = 0; i <= n; ++i)
                    pt[i] = abs(p[i]);
                pt[n] = -pt[n];
                T x = exp((log(-pt[n]) - log(pt[0])) / static_cast<T>(n));
                if (pt[n - 1] != T(0)) {
                    T xm = -pt[n] / pt[n - 1];
                    if (xm < x) x = xm;
                }
                // 区間縮小
                while (true) {
                    T xm = x * T(0.1);
                    T ff = pt[0];
                    for (int i = 1; i <= n; ++i)
                        ff = ff * xm + pt[i];
                    if (ff <= T(0)) break;
                    x = xm;
                }
                T dx = x;
                // Newton 反復で下界を精密化
                while (abs(dx / x) > T(0.005)) {
                    T ff = pt[0], df = ff;
                    for (int i = 1; i < n; ++i) {
                        ff = ff * x + pt[i];
                        df = df * x + ff;
                    }
                    ff = ff * x + pt[n];
                    dx = ff / df;
                    x -= dx;
                }
                T bnd = x;

                // Stage 1: 無シフト K 多項式 (5ステップ)
                int nm1 = n - 1;
                for (int i = 1; i < n; ++i)
                    k[i] = static_cast<T>(n - i) * p[i] / static_cast<T>(n);
                k[0] = p[0];
                T aa = p[n], bb = p[n - 1];
                bool zeroK = (k[n - 1] == T(0));
                for (int jj = 0; jj < 5; ++jj) {
                    T cc = k[n - 1];
                    if (!zeroK) {
                        T t = -aa / cc;
                        for (int i = 0; i < nm1; ++i) {
                            int j = n - i - 1;
                            k[j] = t * k[j - 1] + p[j];
                        }
                        k[0] = p[0];
                        zeroK = (abs(k[n - 1]) <= abs(bb) * eta * T(10));
                    } else {
                        for (int i = 0; i < nm1; ++i) {
                            int j = n - i - 1;
                            k[j] = k[j - 1];
                        }
                        k[0] = T(0);
                        zeroK = (k[n - 1] == T(0));
                    }
                }

                // K 多項式を保存
                for (int i = 0; i < n; ++i)
                    temp[i] = k[i];

                // Stage 2-3: シフト反復 (最大20回のシフト)
                bool found = false;
                for (int cnt = 0; cnt < 20; ++cnt) {
                    T xxx = cosr * xx - sinr * yy;
                    yy = sinr * xx + cosr * yy;
                    xx = xxx;
                    sr = bnd * xx;
                    si = bnd * yy;
                    u = T(-2) * sr;
                    v = bnd;
                    int nz = 0;
                    fxshfr(20 * (cnt + 1), nz);
                    if (nz != 0) {
                        int j = degree - n;
                        zeror[j] = szr;
                        zeroi[j] = szi;
                        foundCount += nz;
                        n -= nz;
                        for (int i = 0; i <= n; ++i)
                            p[i] = qp[i];
                        if (nz == 2) {
                            zeror[j + 1] = lzr;
                            zeroi[j + 1] = lzi;
                        }
                        found = true;
                        break;
                    }
                    // 復元
                    for (int i = 0; i < n; ++i)
                        k[i] = temp[i];
                }
                if (!found) break;  // 20シフトで収束せず
            }

            // 残った1次 or 2次
            if (n == 1) {
                zeror[degree - 1] = -p[1] / p[0];
                zeroi[degree - 1] = T(0);
                foundCount++;
            }
            else if (n == 2) {
                T sr_, si_, lr_, li_;
                quad(p[0], p[1], p[2], sr_, si_, lr_, li_);
                zeror[degree - 2] = sr_;
                zeroi[degree - 2] = si_;
                zeror[degree - 1] = lr_;
                zeroi[degree - 1] = li_;
                foundCount += 2;
            }
        };

        findRoots();
    }

done:
    // 結果を Complex<T> のベクトルに変換
    std::vector<Complex<T>> result;
    result.reserve(foundCount);
    for (int i = 0; i < foundCount; ++i)
        result.push_back(Complex<T>(zeror[i], zeroi[i]));

    return result;
}

// ================================================================
// B-2. Laguerre 法
// ================================================================

/// Laguerre 法による多項式求根
/// 3次収束で1根ずつ求め、デフレーション (因数除去) で次数を下げる
/// 移植元: lib++20 Laguerre.h
template<typename T>
std::vector<Complex<T>> laguerre(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    using std::abs;
    using std::sqrt;

    int n = poly.degree();
    if (n < 1) return {};

    // 複素係数の多項式に変換
    std::vector<Complex<T>> coeffs(n + 1);
    for (int i = 0; i <= n; ++i)
        coeffs[i] = Complex<T>(poly[i]);

    std::vector<Complex<T>> roots;

    // Horner 法で複素多項式を評価
    auto evalPoly = [](const std::vector<Complex<T>>& c, const Complex<T>& z) {
        int deg = static_cast<int>(c.size()) - 1;
        Complex<T> val = c[deg];
        for (int i = deg - 1; i >= 0; --i)
            val = val * z + c[i];
        return val;
    };

    // 微分の評価
    auto evalDeriv = [](const std::vector<Complex<T>>& c, const Complex<T>& z) {
        int deg = static_cast<int>(c.size()) - 1;
        Complex<T> val = c[deg] * T(deg);
        for (int i = deg - 1; i >= 1; --i)
            val = val * z + c[i] * T(i);
        return val;
    };

    // 2階微分の評価
    auto evalDeriv2 = [](const std::vector<Complex<T>>& c, const Complex<T>& z) {
        int deg = static_cast<int>(c.size()) - 1;
        Complex<T> val = c[deg] * T(deg) * T(deg - 1);
        for (int i = deg - 1; i >= 2; --i)
            val = val * z + c[i] * T(i) * T(i - 1);
        return val;
    };

    while (n > 2) {
        T tn = static_cast<T>(n);
        T tn1 = static_cast<T>(n - 1);

        // 初期値: Cauchy の根の上界 r = 1 + max|a_i/a_n| を用いた非自明な点
        T maxRatio = T(0);
        for (int i = 0; i < n; ++i) {
            T ratio = sangi::abs(coeffs[i]) / sangi::abs(coeffs[n]);
            if (ratio > maxRatio) maxRatio = ratio;
        }
        T r0 = T(1) + maxRatio;
        // 実軸上を避けた初期値 (複素数を使うことで複素根も見つかる)
        Complex<T> z(r0 * T(0.4), r0 * T(0.9));

        for (size_t iter = 0; iter < maxIter; ++iter) {
            Complex<T> zo = z;
            Complex<T> pz = evalPoly(coeffs, z);
            Complex<T> dpz = evalDeriv(coeffs, z);
            Complex<T> ddpz = evalDeriv2(coeffs, z);

            // Laguerre の修正公式 (除算回避版)
            // H_tilde = (n-1)^2 * P'(z)^2 - n(n-1) * P(z) * P''(z)
            Complex<T> H = tn1 * (tn1 * (dpz * dpz) - tn * pz * ddpz);
            Complex<T> sqrtH = sangi::sqrt(H);

            Complex<T> d1 = dpz + sqrtH;
            Complex<T> d2 = dpz - sqrtH;

            // 絶対値が大きい方を選択 (数値安定性)
            Complex<T> denom = (normSq(d1) >= normSq(d2)) ? d1 : d2;

            if (normSq(denom) == T(0)) {
                // 分母がゼロ → 摂動を加えてリトライ
                z += Complex<T>(r0 * T(0.1), r0 * T(0.1));
                continue;
            }

            Complex<T> delta = Complex<T>(tn) * pz / denom;
            z -= delta;

            if (sangi::abs(z - zo) < eps)
                break;
        }

        roots.push_back(z);
        n--;

        // デフレーション: P を (x - z) で割る
        std::vector<Complex<T>> newCoeffs(n + 1);
        newCoeffs[n] = coeffs[n + 1];
        for (int i = n - 1; i >= 0; --i)
            newCoeffs[i] = coeffs[i + 1] + newCoeffs[i + 1] * z;
        coeffs = std::move(newCoeffs);
    }

    // 残った1次 or 2次
    if (n == 1) {
        roots.push_back(-coeffs[0] / coeffs[1]);
    }
    else if (n == 2) {
        Complex<T> a = coeffs[2];
        Complex<T> b = coeffs[1];
        Complex<T> c = coeffs[0];
        Complex<T> disc = sangi::sqrt(b * b - Complex<T>(T(4)) * a * c);
        roots.push_back((-b + disc) / (Complex<T>(T(2)) * a));
        roots.push_back((-b - disc) / (Complex<T>(T(2)) * a));
    }

    return roots;
}

// ================================================================
// B-3. Durand-Kerner-Aberth (DKA) 法
// ================================================================

/// Aberth-Ehrlich 法による多項式求根 (全根同時探索, 3次収束)
///
/// 更新式: w[i] = P(z[i]) / P'(z[i])
///         z[i] -= w[i] / (1 - w[i] * Σ_{j≠i} 1/(z[i]-z[j]))
/// 単純な DK 法 (2次収束) と異なり、P'/P の情報を使うため
/// 重根近傍でも収束が安定する。
///
/// 収束後に Newton polish: 元の多項式で各根を精密化する。
template<typename T>
std::vector<Complex<T>> durandKernerAberth(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    using std::abs;
    using std::sqrt;
    using std::cos;
    using std::sin;

    const T pi = std::acos(T(-1));

    int N = poly.degree();
    if (N < 1) return {};
    if (N == 1) return solveLinear(poly);
    if (N == 2) return solveQuadratic(poly);

    // モニック化
    T lc = poly.leadingCoefficient();
    std::vector<T> coeffs(N + 1);
    for (int i = 0; i <= N; ++i)
        coeffs[i] = poly[i] / lc;

    // 根の重心シフト: β = -a_{N-1}/N
    T beta = -coeffs[N - 1] / static_cast<T>(N);

    // シフトした多項式を作る: P(x + β)
    // Horner 法による Taylor シフト
    std::vector<T> shifted(N + 1);
    for (int i = 0; i <= N; ++i)
        shifted[i] = coeffs[i];
    for (int j = 0; j < N; ++j) {
        for (int i = N - 1; i >= j; --i) {
            shifted[i] += beta * shifted[i + 1];
        }
    }

    // 零根を除去
    int zeroRoots = 0;
    while (zeroRoots < N && abs(shifted[zeroRoots]) < eps * T(100)) {
        zeroRoots++;
    }
    int M = N - zeroRoots;  // 有効な次数
    if (M == 0) {
        std::vector<Complex<T>> result(N, Complex<T>(beta));
        return result;
    }

    std::vector<T> workCoeffs(M + 1);
    for (int i = 0; i <= M; ++i)
        workCoeffs[i] = shifted[i + zeroRoots];

    // Horner 法で P(z) と P'(z) を同時評価
    auto evalPP = [&](const Complex<T>& z) -> std::pair<Complex<T>, Complex<T>> {
        Complex<T> p(workCoeffs[M]);
        Complex<T> dp(T(0));
        for (int i = M - 1; i >= 0; --i) {
            dp = dp * z + p;
            p  = p  * z + Complex<T>(workCoeffs[i]);
        }
        return { p, dp };
    };

    // 初期値: Fujiwara 上界半径の円周上に非等間隔配置
    // Fujiwara bound: 2 * max_i( |a_i/a_n|^(1/(n-i)) )
    // workCoeffs はモニック (a_n = 1)
    T r = T(0);
    for (int i = 0; i < M; ++i) {
        T val = std::pow(abs(workCoeffs[i]), T(1) / static_cast<T>(M - i));
        if (val > r) r = val;
    }
    r = std::max(T(1), T(2) * r);

    std::vector<Complex<T>> X(M);
    for (int i = 0; i < M; ++i) {
        // 無理数角度オフセットで初期値の偶然の対称性を回避
        T theta = T(2) * pi * static_cast<T>(i) / static_cast<T>(M)
                  + T(0.4) / static_cast<T>(M);
        X[i] = Complex<T>(r * cos(theta), r * sin(theta));
    }

    // Aberth-Ehrlich 反復
    // 高次多項式では収束に多くの反復を要することがある
    size_t maxAberthIter = std::max(maxIter, static_cast<size_t>(M) * 50);
    T epsSq = eps * eps;
    for (size_t iter = 0; iter < maxAberthIter; ++iter) {
        bool converged = true;
        for (int i = 0; i < M; ++i) {
            auto [pz, dpz] = evalPP(X[i]);

            // Newton 補正: w = P(z) / P'(z)
            T dpNorm = normSq(dpz);
            if (dpNorm < epsSq * epsSq) continue;
            Complex<T> w = pz / dpz;

            // Aberth 補正: Σ 1/(z[i] - z[j])
            Complex<T> aberth_sum(T(0));
            for (int j = 0; j < M; ++j) {
                if (j == i) continue;
                Complex<T> dz = X[i] - X[j];
                T dzNorm = normSq(dz);
                if (dzNorm > epsSq * epsSq)
                    aberth_sum += Complex<T>(T(1)) / dz;
            }

            Complex<T> denom = Complex<T>(T(1)) - w * aberth_sum;
            Complex<T> offset;
            if (normSq(denom) < epsSq * epsSq) {
                // 退化時は Newton ステップにフォールバック
                offset = w;
            } else {
                offset = w / denom;
            }

            // 相対収束判定: |offset| > eps * (|z| + 1)
            T zi_scale = abs(X[i]) + T(1);
            if (abs(offset) > eps * zi_scale)
                converged = false;

            X[i] -= offset;
        }
        if (converged) break;
    }

    // シフトを戻す + 零根を追加
    std::vector<Complex<T>> result;
    result.reserve(N);
    for (int i = 0; i < M; ++i)
        result.push_back(X[i] + Complex<T>(beta));
    for (int i = 0; i < zeroRoots; ++i)
        result.push_back(Complex<T>(beta));

    // Newton polish: 元の多項式で各根を精密化
    Polynomial<T> dp = poly.derivative();
    for (auto& z : result) {
        for (int k = 0; k < 5; ++k) {
            Complex<T> pz = poly(z);
            Complex<T> dpz = dp(z);
            if (normSq(dpz) < epsSq * epsSq) break;
            Complex<T> corr = pz / dpz;
            z -= corr;
            if (normSq(corr) <= epsSq * normSq(z)) break;
        }
        // 実根の虚部を切り落とす
        if (abs(z.im) < eps * (T(1) + abs(z.re)))
            z.im = T(0);
    }

    return result;
}

// ================================================================
// B-4. Bairstow 法
// ================================================================

/// Bairstow 法による多項式求根 (実係数2次因子分解)
///
/// 実係数多項式を反復的に (x² + r·x + s) で合成除算し、
/// Newton 法で r, s を修正して2次因子を抽出する。
/// 実係数のまま計算するため、複素根は共役ペアとして自然に得られる。
///
/// 参考: W.H. Press et al., Numerical Recipes (§9.5)
template<typename T>
std::vector<Complex<T>> bairstow(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    using std::abs;
    using std::sqrt;

    int deg = poly.degree();
    if (deg < 1) return {};
    if (deg == 1) return solveLinear(poly);
    if (deg == 2) return solveQuadratic(poly);

    // 降べきの係数配列 a[0]=最高次, ..., a[n]=定数
    int n = deg;
    std::vector<T> a(n + 1);
    for (int i = 0; i <= n; ++i)
        a[i] = poly[n - i];

    std::vector<Complex<T>> roots;
    roots.reserve(deg);

    while (n > 2) {
        // 初期推定値: 低次係数比に基づく推定
        T r, s;
        if (abs(a[0]) > eps) {
            r = a[1] / a[0];
            s = a[2] / a[0];
        } else {
            r = T(1);
            s = T(1);
        }

        std::vector<T> b(n + 1), c(n + 1);

        for (size_t iter = 0; iter < maxIter; ++iter) {
            // 合成除算: P(x) を (x² + r·x + s) で割る
            // b[j] = a[j] - r·b[j-1] - s·b[j-2]
            b[0] = a[0];
            b[1] = a[1] - r * b[0];
            for (int i = 2; i <= n; ++i)
                b[i] = a[i] - r * b[i - 1] - s * b[i - 2];

            // 偏微分用: 同じ漸化式を b に適用
            // c[j] = b[j] - r·c[j-1] - s·c[j-2]
            // 結果: ∂b[j]/∂r = -c[j-1], ∂b[j]/∂s = -c[j-2]
            c[0] = b[0];
            c[1] = b[1] - r * c[0];
            for (int i = 2; i <= n - 1; ++i)
                c[i] = b[i] - r * c[i - 1] - s * c[i - 2];

            // Newton 系:
            //   [c[n-2]  c[n-3]] [Δr]   [b[n-1]]
            //   [c[n-1]  c[n-2]] [Δs] = [b[n]  ]
            T cn2 = c[n - 2];
            T cn3 = (n >= 3) ? ((n - 3 >= 0) ? c[n - 3] : T(0)) : T(0);
            T cn1 = c[n - 1];

            T det = cn2 * cn2 - cn3 * cn1;
            if (abs(det) < eps * eps) {
                // 特異 → 摂動を加えてやり直す
                r += T(0.5);
                s -= T(0.5);
                continue;
            }

            T dr = (b[n - 1] * cn2 - b[n] * cn3) / det;
            T ds = (b[n] * cn2 - b[n - 1] * cn1) / det;

            r += dr;
            s += ds;

            if (abs(dr) < eps * (T(1) + abs(r)) &&
                abs(ds) < eps * (T(1) + abs(s))) {
                break;
            }
        }

        // 2次因子 x² + r·x + s の根を求める
        T disc = r * r - T(4) * s;
        if (disc >= T(0)) {
            T sq = sqrt(disc);
            roots.push_back(Complex<T>((-r + sq) / T(2)));
            roots.push_back(Complex<T>((-r - sq) / T(2)));
        } else {
            T sq = sqrt(-disc);
            roots.push_back(Complex<T>(-r / T(2), sq / T(2)));
            roots.push_back(Complex<T>(-r / T(2), -sq / T(2)));
        }

        // デフレーション: 商多項式 b[0..n-2] を次の a にする
        n -= 2;
        std::vector<T> newA(n + 1);
        for (int i = 0; i <= n; ++i)
            newA[i] = b[i];
        a = std::move(newA);
    }

    // 残った1次 or 2次
    if (n == 1) {
        roots.push_back(Complex<T>(-a[1] / a[0]));
    } else if (n == 2) {
        T disc = a[1] * a[1] - T(4) * a[0] * a[2];
        if (disc >= T(0)) {
            T sq = sqrt(disc);
            roots.push_back(Complex<T>((-a[1] + sq) / (T(2) * a[0])));
            roots.push_back(Complex<T>((-a[1] - sq) / (T(2) * a[0])));
        } else {
            T sq = sqrt(-disc);
            roots.push_back(Complex<T>(-a[1] / (T(2) * a[0]), sq / (T(2) * a[0])));
            roots.push_back(Complex<T>(-a[1] / (T(2) * a[0]), -sq / (T(2) * a[0])));
        }
    }

    return roots;
}

// ================================================================
// C. 統合ディスパッチャ
// ================================================================

/// 多項式の全根を求める (次数に応じて最適なアルゴリズムを自動選択)
///
/// アルゴリズム選択:
///   1-4次: 公式解 (桁落ち対策付き)
///   5-19次: Jenkins-Traub (1根ずつ+デフレーション, 高精度)
///   20次以上: Aberth-Ehrlich + Newton polish (全根同時探索,
///             デフレーション誤差蓄積を回避)
template<typename T>
std::vector<Complex<T>> solvePolynomial(
    const Polynomial<T>& p,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    // 先頭の零係数を除去して実効次数を求める
    int deg = p.degree();
    if (deg < 0) return {};  // 零多項式

    // 実効的な多項式を取得 (正規化済みなので degree() が正しい)
    switch (deg) {
        case 0: return {};  // 定数多項式 → 根なし
        case 1: return solveLinear(p);
        case 2: return solveQuadratic(p);
        case 3: return solveCubic(p);
        case 4: return solveQuartic(p);
        default:
            // 20次以上: DKA (全根同時 + Newton polish)
            // デフレーション誤差蓄積を回避
            if (deg >= 20)
                return durandKernerAberth(p, eps, maxIter);
            // 5-19次: Jenkins-Traub (1根ずつ、高精度)
            return jenkinsTraub(p);
    }
}

// ================================================================
// D. Vincent 法による実根分離 (Vincent Real Root Isolation)
// ================================================================

/**
 * @brief 実根分離区間
 */
template<typename T>
struct RealRootInterval {
    T lower;           ///< 区間の下限
    T upper;           ///< 区間の上限
    bool exact;        ///< true なら lower == upper (厳密な根)
};

namespace detail {

/// 係数列の符号変動数 (Descartes の符号規則)
template<typename T>
std::size_t signVariations(const std::vector<T>& coeffs) {
    std::size_t count = 0;
    int last_sign = 0;  // -1, 0, +1
    for (const auto& c : coeffs) {
        int s = (c > T(0)) ? 1 : (c < T(0)) ? -1 : 0;
        if (s == 0) continue;
        if (last_sign != 0 && s != last_sign) ++count;
        last_sign = s;
    }
    return count;
}

/// 多項式の符号変動数
template<typename T>
std::size_t signVariations(const Polynomial<T>& p) {
    return signVariations(p.coefficients());
}

/// Taylor シフト: p(x + c) を計算 (Horner 方式, O(n²))
template<typename T>
Polynomial<T> taylorShift(const Polynomial<T>& p, const T& c) {
    int n = p.degree();
    if (n < 0) return p;

    // coeffs[i] を p(x+c) の係数に変換
    std::vector<T> a(p.coefficients());
    for (int i = 0; i < n; ++i) {
        for (int j = n - 1; j >= i; --j) {
            a[j] += c * a[j + 1];
        }
    }
    return Polynomial<T>(std::move(a));
}

/// 逆数多項式: x^n * p(1/x) (係数を逆順にする)
template<typename T>
Polynomial<T> reciprocalPoly(const Polynomial<T>& p) {
    auto c = p.coefficients();
    std::reverse(c.begin(), c.end());
    return Polynomial<T>(std::move(c));
}

/// Cauchy 上界: 正の実根の上界 max(1, Σ|a_i/a_n|)
template<typename T>
T cauchyBound(const Polynomial<T>& p) {
    int n = p.degree();
    if (n <= 0) return T(1);
    T lc = std::abs(p[n]);
    T bound = T(0);
    for (int i = 0; i < n; ++i)
        bound = std::max(bound, std::abs(p[i]) / lc);
    return T(1) + bound;
}

/// 無平方化: p / gcd(p, p')
template<typename T>
Polynomial<T> squareFree(const Polynomial<T>& p) {
    auto dp = p.derivative();
    if (dp.isZero()) return p;
    auto g = gcd(p, dp);
    if (g.degree() <= 0) return p;
    return p / g;
}

/// 二分法で根を精密化
template<typename T>
T bisectRoot(const Polynomial<T>& p, T lo, T hi, T tol, std::size_t max_iter = 100) {
    T flo = p(lo), fhi = p(hi);
    if (std::abs(flo) <= tol) return lo;
    if (std::abs(fhi) <= tol) return hi;
    // 符号が同じなら中点を返す
    if ((flo > T(0)) == (fhi > T(0))) return (lo + hi) / T(2);
    for (std::size_t i = 0; i < max_iter; ++i) {
        T mid = (lo + hi) / T(2);
        if (hi - lo < tol) return mid;
        T fm = p(mid);
        if (std::abs(fm) <= tol) return mid;
        if ((fm > T(0)) == (flo > T(0))) { lo = mid; flo = fm; }
        else { hi = mid; fhi = fm; }
    }
    return (lo + hi) / T(2);
}

/// 正の実根の分離 (VCA bisection)
/// Möbius 変換 x = (a*t + b)/(c*t + d), t ∈ (0, ∞) を追跡
template<typename T>
void vincentPositiveRoots(
    const Polynomial<T>& p_orig,
    T upper_bound,
    std::vector<RealRootInterval<T>>& result,
    T epsilon)
{
    struct WorkItem {
        Polynomial<T> q;  // 変換後の多項式
        T a, b, c, d;     // Möbius パラメータ
    };

    std::vector<WorkItem> stack;
    stack.push_back({p_orig, T(1), T(0), T(0), T(1)});

    std::size_t max_work = static_cast<std::size_t>(p_orig.degree()) * 500;
    std::size_t work = 0;

    while (!stack.empty() && work < max_work) {
        ++work;
        auto [q, a, b, c, d] = std::move(stack.back());
        stack.pop_back();

        int deg = q.degree();
        if (deg < 1) continue;

        // 主係数を正に正規化
        if (q[deg] < T(0)) q = q * T(-1);

        std::size_t v = signVariations(q);
        if (v == 0) continue;

        if (v == 1) {
            // q の正の実根を求め、Möbius 逆変換で元の座標に戻す
            T t_root;
            if (deg == 1) {
                // 1 次: 直接解く
                t_root = -q[0] / q[1];
            } else {
                // bisection で q の正の根を求める
                T ub = cauchyBound(q);
                t_root = bisectRoot(q, T(0), ub,
                                    std::numeric_limits<T>::epsilon() * T(100));
            }
            if (t_root >= T(0)) {
                T x_root = (a * t_root + b) / (c * t_root + d);
                T margin = std::max(std::abs(x_root) * epsilon, epsilon);
                result.push_back({x_root - margin, x_root + margin, false});
            }
            continue;
        }

        // t = 1 での根をチェック (= x = (a+b)/(c+d))
        T val_at_1 = q(T(1));
        if (std::abs(val_at_1) <= epsilon * (std::abs(q[deg]) + T(1))) {
            T exact_root = (a + b) / (c + d);
            result.push_back({exact_root, exact_root, true});
            // q / (x - 1) でデフレーション
            auto [quot, rem] = q.divmod(Polynomial<T>({T(-1), T(1)}));
            q = std::move(quot);
            if (q.degree() < 1) continue;
            if (q[q.degree()] < T(0)) q = q * T(-1);
            v = signVariations(q);
            if (v == 0) continue;
        }

        // (1, ∞) 部分: t → t + 1, q₁(t) = q(t + 1)
        // Möbius: (a(t+1)+b)/(c(t+1)+d) = (at+(a+b))/(ct+(c+d))
        auto q1 = taylorShift(q, T(1));
        stack.push_back({std::move(q1), a, a + b, c, c + d});

        // (0, 1) 部分: t → 1/(t+1)
        // q₂(t) = (t+1)^n * q(1/(t+1)) = taylorShift(reciprocal(q), 1)
        // Möbius: (bt+(a+b))/(dt+(c+d))
        auto q2 = taylorShift(reciprocalPoly(q), T(1));
        stack.push_back({std::move(q2), b, a + b, d, c + d});
    }
}

} // namespace detail

/**
 * @brief Vincent 法による実根分離
 *
 * 多項式の全実根を互いに素な区間に分離する。
 * Vincent-Collins-Akritas (VCA) アルゴリズムの bisection 版を使用。
 *
 * 各区間にはちょうど 1 つの実根が含まれる。
 * exact == true の場合は厳密な根値 (lower == upper)。
 *
 * @param p       入力多項式
 * @param epsilon 区間幅の許容下限 (デフォルト: 100ε)
 * @return 分離区間のベクトル (lower ≤ upper, ソート済み)
 */
template<typename T>
[[nodiscard]]
std::vector<RealRootInterval<T>> vincentRealRootIsolation(
    const Polynomial<T>& p,
    T epsilon = std::numeric_limits<T>::epsilon() * T(100))
{
    int deg = p.degree();
    if (deg <= 0) return {};

    std::vector<RealRootInterval<T>> result;

    // 1次は直接解く
    if (deg == 1) {
        T root = -p[0] / p[1];
        result.push_back({root, root, true});
        return result;
    }

    // x = 0 の根を抽出
    Polynomial<T> q = p;
    int zero_roots = 0;
    while (q.degree() >= 1 && std::abs(q[0]) <= epsilon * std::abs(q[q.degree()])) {
        // q / x
        std::vector<T> c(q.coefficients().begin() + 1, q.coefficients().end());
        q = Polynomial<T>(std::move(c));
        ++zero_roots;
    }
    if (zero_roots > 0) {
        result.push_back({T(0), T(0), true});
    }

    if (q.degree() < 1) {
        return result;
    }

    // 無平方化
    auto sf = detail::squareFree(q);
    T bound = detail::cauchyBound(sf);

    // 正の根を分離
    detail::vincentPositiveRoots(sf, bound, result, epsilon);

    // 負の根: p(-x) の正の根を分離し、符号を反転
    // p(-x): 奇数次係数の符号を反転
    auto coeffs = sf.coefficients();
    for (std::size_t i = 0; i < coeffs.size(); ++i) {
        if (i % 2 == 1) coeffs[i] = -coeffs[i];
    }
    Polynomial<T> neg_poly(std::move(coeffs));
    T neg_bound = detail::cauchyBound(neg_poly);

    std::vector<RealRootInterval<T>> neg_intervals;
    detail::vincentPositiveRoots(neg_poly, neg_bound, neg_intervals, epsilon);
    for (auto& iv : neg_intervals) {
        result.push_back({-iv.upper, -iv.lower, iv.exact});
    }

    // ソート
    std::sort(result.begin(), result.end(),
        [](const RealRootInterval<T>& a, const RealRootInterval<T>& b) {
            return a.lower < b.lower;
        });

    return result;
}

// ================================================================
// Bernoulli 法 (冪乗法による支配的根の抽出)
// ================================================================

/**
 * @brief Bernoulli 法 (冪乗法) で多項式の絶対値最大の根 (支配的根) を求める
 *
 * 多項式 p(x) = a_n x^n + ... + a_1 x + a_0 をモニック化し、
 * 特性方程式に基づく線形漸化式 y_k = -c_{n-1} y_{k-1} - ... - c_0 y_{k-n}
 * のインパルス応答の比 y_k / y_{k-1} が支配的根に収束する性質を利用する。
 *
 * 参考: 戸川隼人「数値計算」p74
 *
 * @param poly 多項式 (次数 ≥ 1)
 * @param tolerance 連続する近似根の差がこれ以下で収束
 * @param max_iterations 最大反復回数
 * @return 支配的根と収束状態
 *
 * @note 絶対値最大の根が重根または複素共役対で等絶対値の場合は収束しない。
 */
template<typename T>
std::pair<T, bool> bernoulli_method(
    const Polynomial<T>& poly,
    T tolerance = std::numeric_limits<T>::epsilon() * T(1000),
    size_t max_iterations = 1000) {

    int n = poly.degree();
    if (n <= 0) return {T(0), false};

    // モニック化: 最高次係数で割る
    T lead = poly[n];
    std::vector<T> c(n);
    for (int i = 0; i < n; ++i) {
        c[i] = poly[i] / lead;
    }

    // 状態ベクトル y[0..n-1]: インパルス応答の最新 n 個
    // 初期値: y[0]=1 (インパルス), y[1..n-1]=0
    std::vector<T> y(n, T(0));
    y[0] = T(1);

    // 線形漸化式でインパルス応答を計算
    // y_new = -(c[n-1]*y[0] + c[n-2]*y[1] + ... + c[0]*y[n-1])
    T r_prev = T(0);
    bool converged = false;
    T r = T(0);

    // 最初の数ステップでは y[previous] = 0 になりうるのでスキップ
    size_t warmup = static_cast<size_t>(n) + 2;

    for (size_t iter = 0; iter < max_iterations + warmup; ++iter) {
        T y_new = T(0);
        for (int j = 0; j < n; ++j) {
            y_new -= c[n - 1 - j] * y[j];
        }

        // 比率計算 (ウォームアップ後)
        if (iter >= warmup && std::abs(y[0]) > std::numeric_limits<T>::epsilon()) {
            r = y_new / y[0];

            if (iter > warmup && std::abs(r - r_prev) <= tolerance) {
                converged = true;
                break;
            }
            r_prev = r;
        }

        // シフト: y[n-1] を捨てて y_new を先頭に
        for (int j = n - 1; j > 0; --j) {
            y[j] = y[j - 1];
        }
        y[0] = y_new;

        // オーバーフロー防止: 値が大きくなりすぎたらスケーリング
        T max_val = T(0);
        for (int j = 0; j < n; ++j) {
            T av = std::abs(y[j]);
            if (av > max_val) max_val = av;
        }
        if (max_val > T(1e100)) {
            for (int j = 0; j < n; ++j) {
                y[j] /= max_val;
            }
        }
    }

    return {r, converged};
}

} // namespace sangi

#endif // SANGI_POLYNOMIAL_ROOTS_HPP