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

// polynomial_roots.hpp
// Numerical root-finding algorithms for polynomials
//
// Low degree (1 to 4): closed-form solutions robust against cancellation
// High degree (5 and above): Jenkins-Traub, Laguerre, DKA iterative methods
//
// All functions return a vector of Complex<T> (real roots have 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 {

// ================================================================
// Internal helpers
// ================================================================

namespace detail {

// Sign function
template<typename T>
T sgn(T val) {
    if (val > T(0)) return T(1);
    if (val < T(0)) return T(-1);
    return T(0);
}

// Cube root (real)
template<typename T>
T cubeRoot(T x) {
    using std::cbrt;
    return cbrt(x);
}

} // namespace detail

// ================================================================
// A-1. solveLinear (linear equation)
// ================================================================

/// Returns the solution of the linear equation a*x + b = 0
/// p.coefficients(): [b, a] (ascending powers)
template<typename T>
std::vector<Complex<T>> solveLinear(const Polynomial<T>& p) {
    assert(p.degree() == 1);

    T a = p[1];  // coefficient of x
    T b = p[0];  // constant term

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

// ================================================================
// A-2. solveQuadratic (quadratic equation)
// ================================================================

/// Returns the solution of the quadratic equation a*x^2 + b*x + c = 0
/// Cancellation mitigation: compute the larger-magnitude root first, then the smaller root via Vieta's formulas
/// [Source] Okumura, "Encyclopedia of Algorithms", p205
template<typename T>
std::vector<Complex<T>> solveQuadratic(const Polynomial<T>& p) {
    assert(p.degree() == 2);

    T a = p[2];  // coefficient of x^2
    T b = p[1];  // coefficient of x
    T c = p[0];  // constant term

    // Make monic: 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;  // discriminant D/4

    if (disc > T(0)) {
        // 2 real roots — cancellation mitigation
        // Compute the larger-magnitude root first
        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)) {
        // Complex conjugate roots
        T re = -halfP;
        T im = sqrt(-disc);
        return { Complex<T>(re, im), Complex<T>(re, -im) };
    }
    else {
        // Double root
        T x = -halfP;
        return { Complex<T>(x), Complex<T>(x) };
    }
}

// ================================================================
// A-3. solveCubic (cubic equation)
// ================================================================

/// Returns the solution of the cubic equation a*x^3 + b*x^2 + c*x + d = 0
/// Cardano's formula + trigonometric method (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];

    // Make monic
    b /= a;  c /= a;  d /= a;

    // Tschirnhaus transformation: 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

    // Discriminant Δ = -(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 real roots → trigonometric method
        T r = sqrt(-pp / T(3));                        // r = √(-p/3)
        T cosTheta = -qq / (T(2) * r * r * r);        // cos(θ) = -q/(2r^3)
        // Clamp when numerical error pushes it beyond ±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 real root + complex conjugate roots
        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));

        // Remaining 2 roots: divide by (x - realRoot) and solve the quadratic
        // From t^3 + pt + q = (t - (A+B))(t^2 + (A+B)t + ...),
        // obtain the quadratic coefficients via synthetic division
        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: repeated root(s)
        if (pp == T(0) && qq == T(0)) {
            // Triple root: 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 {
            // Double root + single root
            // When t³+pt+q=0 is (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 (quartic equation)
// ================================================================

/// Returns the solution of the quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
/// Ferrari's method: factor into two quadratics using a real root of the resolvent cubic
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];

    // Make monic
    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
        // Letting 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;
    }

    // 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);

    // Select the largest among the real roots (numerical stability)
    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;
        }
    }

    // Factor into two quadratics
    // 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 {
        // Fallback: tiny negative value due to numerical error
        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 {
        // Case α ≈ 0
        T sqrtY0sq_r = sqrt(abs(y0 * y0 - rr));
        beta  = y0 - sqrtY0sq_r;
        gamma = y0 + sqrtY0sq_r;
    }

    // Solve the two quadratic equations
    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);

    // Undo the shift
    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 method
// ================================================================

/// Polynomial root-finding by the Jenkins-Traub method (based on TOMS493)
/// The most robust polynomial root-finding method. A three-stage algorithm.
/// Based on C. Bond's implementation
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 {};

    // Build the coefficient array in descending order (TOMS493 convention)
    // poly is in ascending powers: coeffs_[0]=constant, coeffs_[n]=leading coefficient
    // op[0]=leading coefficient, op[n]=constant term
    std::vector<T> op(degree + 1);
    for (int i = 0; i <= degree; ++i)
        op[i] = poly[degree - i];

    // Invalid if the leading coefficient is zero
    if (op[0] == T(0))
        return {};

    // Storage for results
    std::vector<T> zeror(degree), zeroi(degree);
    int foundCount = 0;

    // Working arrays
    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);

    // Machine constants
    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;

    // Rotation angle (94 degrees)
    T rot = T(94) * T(0.017453293);
    T cosr = cos(rot);
    T sinr = sin(rot);
    T xx = sqrt(T(0.5));
    T yy = -xx;

    // Remove zero roots first
    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;

    // Copy the coefficients into p
    for (int i = 0; i <= n; ++i)
        p[i] = op[i];

    // Main root-finding loop
    {
        // Local variables
        T sr, si, u, v, a, b, c, d;
        T a1, a3, a7, ee, f, g, h;
        T szr, szi, lzr, lzi;

        // --- Define the internal functions as lambdas ---

        // quadsd: synthetic division of p by the quadratic (1, u, v)
        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: solution of the quadratic equation 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;
            }
            // Avoid overflow in the discriminant
            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)) {
                // Complex conjugate roots
                sr_out = -bb / qa;
                lr_out = sr_out;
                si_out = abs(dd / qa);
                li_out = -si_out;
            } else {
                // Real roots
                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: compute the scalar quantities
        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: compute the next K polynomial
        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: estimate the new 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: iteration for a real zero
        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;
                // Update the K polynomial
                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: iteration for a quadratic factor
        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);
                // Upper bound on the round-off error
                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: fixed-shift iteration
        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;

                        // Try whichever converges faster
                        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;
            }
        };

        // --- Main loop ---
        auto findRoots = [&]() {
            while (n > 2) {
                // Coefficient scaling
                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;
                    }
                }

                // Compute a lower bound on the roots
                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;
                }
                // Shrink the interval
                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;
                // Refine the lower bound by Newton iteration
                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: no-shift K polynomial (5 steps)
                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));
                    }
                }

                // Save the K polynomial
                for (int i = 0; i < n; ++i)
                    temp[i] = k[i];

                // Stage 2-3: shifted iteration (up to 20 shifts)
                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;
                    }
                    // Restore
                    for (int i = 0; i < n; ++i)
                        k[i] = temp[i];
                }
                if (!found) break;  // did not converge within 20 shifts
            }

            // Remaining linear or quadratic
            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:
    // Convert the result into a vector of 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 method
// ================================================================

/// Polynomial root-finding by the Laguerre method
/// Finds one root at a time with cubic convergence, lowering the degree by deflation (factor removal)
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 {};

    // Convert to a complex-coefficient polynomial
    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;

    // Evaluate the complex polynomial by Horner's method
    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;
    };

    // Evaluate the derivative
    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;
    };

    // Evaluate the second derivative
    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);

        // Initial value: a nontrivial point using the Cauchy root bound 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;
        // Initial value off the real axis (using a complex number also lets complex roots be found)
        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 correction formula (division-avoiding form)
            // 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;

            // Choose the larger-magnitude one (numerical stability)
            Complex<T> denom = (normSq(d1) >= normSq(d2)) ? d1 : d2;

            if (normSq(denom) == T(0)) {
                // Denominator is zero → add a perturbation and retry
                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--;

        // Deflation: divide P by (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);
    }

    // Remaining linear or quadratic
    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) method
// ================================================================

/// Polynomial root-finding by the Aberth-Ehrlich method (all roots simultaneously, cubic convergence)
///
/// Update formula: w[i] = P(z[i]) / P'(z[i])
///         z[i] -= w[i] / (1 - w[i] * Σ_{j≠i} 1/(z[i]-z[j]))
/// Unlike the plain DK method (quadratic convergence), it uses P'/P information,
/// so convergence stays stable even near repeated roots.
///
/// After convergence, Newton polish: refine each root with the original polynomial.
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);

    // Make monic
    T lc = poly.leadingCoefficient();
    std::vector<T> coeffs(N + 1);
    for (int i = 0; i <= N; ++i)
        coeffs[i] = poly[i] / lc;

    // Centroid shift of the roots: β = -a_{N-1}/N
    T beta = -coeffs[N - 1] / static_cast<T>(N);

    // Build the shifted polynomial: P(x + β)
    // Taylor shift by Horner's method
    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];
        }
    }

    // Remove zero roots
    int zeroRoots = 0;
    while (zeroRoots < N && abs(shifted[zeroRoots]) < eps * T(100)) {
        zeroRoots++;
    }
    int M = N - zeroRoots;  // effective degree
    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];

    // Evaluate P(z) and P'(z) simultaneously by Horner's method
    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 };
    };

    // Initial values: placed non-uniformly on the circle of the Fujiwara upper-bound radius
    // Fujiwara bound: 2 * max_i( |a_i/a_n|^(1/(n-i)) )
    // workCoeffs is monic (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) {
        // Avoid accidental symmetry in the initial values via an irrational angle offset
        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 iteration
    // High-degree polynomials may require many iterations to converge
    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 correction: w = P(z) / P'(z)
            T dpNorm = normSq(dpz);
            if (dpNorm < epsSq * epsSq) continue;
            Complex<T> w = pz / dpz;

            // Aberth correction: Σ 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) {
                // On degeneracy, fall back to the Newton step
                offset = w;
            } else {
                offset = w / denom;
            }

            // Relative convergence test: |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;
    }

    // Undo the shift + add the zero roots
    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: refine each root with the original polynomial
    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;
        }
        // Drop the imaginary part of a real root
        if (abs(z.im) < eps * (T(1) + abs(z.re)))
            z.im = T(0);
    }

    return result;
}

// ================================================================
// B-3b. Durand-Kerner (Weierstrass) method
// ================================================================

/// Polynomial root-finding by the Durand-Kerner (Weierstrass) method (all roots simultaneously, classical)
///
/// Updates each root simultaneously by z_i ← z_i − P(z_i) / (lc·Π_{j≠i}(z_i − z_j)).
/// A classical simultaneous iteration without the Aberth correction, more naive than durandKernerAberth
/// (quadratic convergence). Starting from the complex initial values (0.4+0.9i)^i, it converges to all roots when there are no repeated roots.
///
/// References: Durand (1960), Kerner (1966), Weierstrass (1891)
template<typename T>
std::vector<Complex<T>> weierstrass(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000)
{
    using std::abs;
    int N = poly.degree();
    if (N < 1) return {};
    if (N == 1) return solveLinear(poly);
    if (N == 2) return solveQuadratic(poly);

    const T lc = poly.leadingCoefficient();

    // Initial values: (0.4 + 0.9i)^i  (a classical seed: non-real, off the unit circle)
    std::vector<Complex<T>> z(N);
    Complex<T> seed(T(0.4), T(0.9));
    Complex<T> cur(T(1));
    for (int i = 0; i < N; ++i) { cur = cur * seed; z[i] = cur; }

    for (size_t iter = 0; iter < maxIter; ++iter) {
        T maxCorr = T(0);
        for (int i = 0; i < N; ++i) {       // Gauss-Seidel (apply updates immediately)
            Complex<T> denom(T(1));
            for (int j = 0; j < N; ++j)
                if (j != i) denom = denom * (z[i] - z[j]);
            if (abs(denom) < eps * eps) continue;   // roots too close together
            Complex<T> w = poly(z[i]) / (Complex<T>(lc) * denom);   // includes making monic
            z[i] = z[i] - w;
            T c = abs(w);
            if (c > maxCorr) maxCorr = c;
        }
        if (maxCorr < eps) break;
    }

    // Drop the imaginary part of real roots
    for (auto& zi : z)
        if (abs(zi.im) < eps * (T(1) + abs(zi.re))) zi.im = T(0);
    return z;
}

// ================================================================
// B-3c. Graeffe root-squaring method
// ================================================================

/// Polynomial root-finding by Graeffe's root-squaring method (for real roots with well-separated moduli)
///
/// Repeatedly forms p(x)·p(−x) = s(x²), squaring the roots r_i → r_i^{2^m}.
/// For roots with separated moduli, the monic coefficient ratio after squaring collapses to the dominant term, so
/// the modulus can be read off as |r_i| ≈ |a_{n−i}/a_{n−i+1}|^{1/2^m}. The sign of a real root is
/// recovered from P(±|r_i|), and finally polished with a few Newton steps.
///
/// Note (scope): equal-modulus roots, complex conjugate pairs, and repeated roots cannot be separated by
/// basic Graeffe (this implementation targets real roots with separated moduli). Because coefficients grow
/// doubly exponentially, each squaring is monic-normalized, and the number of squarings is limited by maxSquarings. A historical method.
///
/// @param maxSquarings number of squarings (default 6 = exponent 2^6=64)
template<typename T>
std::vector<Complex<T>> graeffe(
    const Polynomial<T>& poly,
    int maxSquarings = 6,
    T eps = std::numeric_limits<T>::epsilon() * T(100))
{
    using std::abs;
    int n = poly.degree();
    if (n < 1) return {};
    if (n == 1) return solveLinear(poly);
    if (n == 2) return solveQuadratic(poly);

    const T tiny = std::numeric_limits<T>::epsilon() * T(10);

    // Monic coefficients a[0..n], a[n]=1
    std::vector<T> a(n + 1);
    const T lc = poly.leadingCoefficient();
    for (int i = 0; i <= n; ++i) a[i] = poly[i] / lc;

    int sq = 0;
    for (sq = 0; sq < maxSquarings; ++sq) {
        // r(x) = a(x)·a(−x) = Σ_{i,j} a_i a_j (−1)^j x^{i+j}
        std::vector<T> r(2 * n + 1, T(0));
        for (int i = 0; i <= n; ++i) {
            T ai = a[i];
            if (ai == T(0)) continue;
            for (int j = 0; j <= n; ++j) {
                T sign = (j % 2 == 0) ? T(1) : T(-1);
                r[i + j] += ai * a[j] * sign;
            }
        }
        // s(y) = even-degree extraction of r (y = x²)
        std::vector<T> b(n + 1);
        for (int k = 0; k <= n; ++k) b[k] = r[2 * k];
        T bn = b[n];
        if (std::abs(bn) < tiny) break;             // degenerate
        for (int k = 0; k <= n; ++k) b[k] /= bn;    // make monic (b[n]=±1)

        // Overflow guard: coefficients grow doubly exponentially, so stop before exceeding the double range.
        // Making monic normalizes only the leading coefficient while the others grow, so squaring further would yield inf.
        T mx = T(0);
        for (int k = 0; k <= n; ++k) mx = std::max(mx, std::abs(b[k]));
        if (!std::isfinite(mx) || mx > T(1e150)) break;   // keep a and sq at the previous stage

        a = b;
    }

    // |r_i| = |a_{n-i}/a_{n-i+1}|^{1/2^sq}
    const T E = std::pow(T(2), T(sq));
    std::vector<Complex<T>> roots;
    roots.reserve(n);
    for (int i = 1; i <= n; ++i) {
        T ratio = std::abs(a[n - i]) / (std::abs(a[n - i + 1]) + tiny);
        T mag = std::pow(ratio, T(1) / E);

        // Sign recovery: whichever of P(+mag) and P(-mag) is smaller
        T r0 = (std::abs(poly(mag)) <= std::abs(poly(-mag))) ? mag : -mag;

        // Polish with a few Newton steps (on the original polynomial)
        Polynomial<T> dp = poly.derivative();
        for (int k = 0; k < 8; ++k) {
            T pv = poly(r0), dpv = dp(r0);
            if (std::abs(dpv) < tiny) break;
            T step = pv / dpv;
            r0 -= step;
            if (std::abs(step) < eps * (T(1) + std::abs(r0))) break;
        }
        roots.push_back(Complex<T>(r0));
    }
    return roots;
}

// ================================================================
// B-3d. Lehmer-Schur method
// ================================================================

namespace detail_lehmer {

// Whether there is a root inside the unit disk |w|<1 (Lehmer's Schur-Cohn test, boolean).
//   q: complex coefficients (ascending powers). Iterate the Schur transform Tq = conj(a0)·q − an·q*
//   until the degree drops; if |a0|²−|an|²<0 occurs along the way, there is a root inside. Degeneracy (|a0|≈|an|) is
//   conservatively treated as true (so as not to miss roots).
template<typename T>
bool hasRootInUnitDisk(std::vector<Complex<T>> q, T tol) {
    using std::abs;
    while (q.size() >= 2) {
        while (q.size() >= 2 && normSq(q.back()) <= tol * tol) q.pop_back();
        if (q.size() < 2) break;
        const std::size_t n = q.size() - 1;
        T scale = T(0);
        for (auto& c : q) scale = std::max(scale, normSq(c));
        if (scale <= tol * tol) break;
        T d = normSq(q[0]) - normSq(q[n]);     // |a0|² − |an|²
        if (d < -tol * scale) return true;       // root inside
        if (d <= tol * scale) return true;       // degenerate → conservatively true
        // Schur transform T q = conj(a0)·q − an·q*  (degree n−1)
        std::vector<Complex<T>> tq(n);
        const Complex<T> a0 = q[0], an = q[n];
        for (std::size_t k = 0; k < n; ++k)
            tq[k] = conj(a0) * q[k] - an * conj(q[n - k]);
        q.swap(tq);
    }
    return false;
}

// Returns the complex coefficients (ascending powers) of P(c + r·w) via Taylor shift + scale.
template<typename T>
std::vector<Complex<T>> shiftScale(const std::vector<Complex<T>>& a,
                                   Complex<T> c, T r) {
    const std::size_t n = a.size() - 1;
    std::vector<Complex<T>> b = a;
    for (std::size_t i = 0; i < n; ++i)                 // Taylor shift P(c+t)
        for (std::size_t j = n - 1; ; --j) {
            b[j] = b[j] + c * b[j + 1];
            if (j == i) break;
        }
    Complex<T> rk(T(1));                                 // scale by t = r·w
    for (std::size_t k = 0; k <= n; ++k) { b[k] = b[k] * rk; rk = rk * Complex<T>(r); }
    return b;
}

// Evaluate a complex-coefficient polynomial by Horner.
template<typename T>
Complex<T> evalComplex(const std::vector<Complex<T>>& a, Complex<T> z) {
    Complex<T> p = a.back();
    for (std::size_t i = a.size() - 1; i-- > 0; ) p = p * z + a[i];
    return p;
}

// Newton polish on a complex-coefficient polynomial.
template<typename T>
Complex<T> newtonPolish(const std::vector<Complex<T>>& a, Complex<T> z, int iters, T tol) {
    const std::size_t n = a.size() - 1;
    for (int it = 0; it < iters; ++it) {
        Complex<T> p = a[n], dp(T(0));
        for (std::size_t i = n; i-- > 0; ) { dp = dp * z + p; p = p * z + a[i]; }
        if (normSq(dp) <= tol * tol) break;
        Complex<T> step = p / dp;
        z = z - step;
        if (normSq(step) <= tol * tol * (T(1) + normSq(z))) break;
    }
    return z;
}

// Synthetic division of a complex-coefficient polynomial by (w − root) (deflation).
template<typename T>
std::vector<Complex<T>> deflate(const std::vector<Complex<T>>& a, Complex<T> root) {
    const std::size_t n = a.size() - 1;
    std::vector<Complex<T>> b(n);
    b[n - 1] = a[n];
    for (std::size_t k = n - 1; k-- > 0; ) b[k] = a[k + 1] + root * b[k + 1];
    return b;
}

}  // namespace detail_lehmer

/// Polynomial root-finding by the Lehmer-Schur method (disk subdivision in the complex plane)
///
/// Starting from a disk containing all roots (Cauchy bound), cover the disk with 9 sub-disks of half the radius
/// (center + 8 surrounding), apply the Schur-Cohn test (Lehmer's unit-disk
/// test) to each sub-disk, and descend into the sub-disk that contains a root. Once the radius is small enough, use the center
/// as a seed for Newton polish → deflation, repeating this degree-many times to obtain all roots.
///
/// The test is based on the exact Schur-Cohn transform. Since seed accuracy comes from the subdivision and root accuracy from Newton polish,
/// even if the cover has some gaps, deflation can still extract all roots.
///
/// References: Lehmer (1961), Schur-Cohn test
template<typename T>
std::vector<Complex<T>> lehmerSchur(
    const Polynomial<T>& poly,
    T tol = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxDescend = 200)
{
    using std::abs;
    using namespace detail_lehmer;

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

    std::vector<Complex<T>> orig(N + 1), work;
    for (int i = 0; i <= N; ++i) orig[i] = Complex<T>(poly[i]);
    work = orig;

    // Cauchy bound: all roots lie in |z| < R
    T an = abs(orig[N]);
    if (!(an > T(0))) return {};
    T m = T(0);
    for (int i = 0; i < N; ++i) m = std::max(m, abs(orig[i]));
    const T R = T(1) + m / an;
    const T seedTol = std::max(tol, T(1e-8) * R);

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

    for (int found = 0; found < N; ++found) {
        const std::size_t deg = work.size() - 1;
        if (deg == 1) {                                  // linear: solve directly
            roots.push_back(Complex<T>(T(0)) - work[0] / work[1]);
            break;
        }

        // Descend by disk subdivision into the sub-disk containing a root
        Complex<T> c(T(0));
        T r = R;
        for (size_t step = 0; step < maxDescend && r > seedTol; ++step) {
            // Cover: center (c, r/2) + 8 surrounding (c + 0.75r·e^{iπk/4}, r/2)
            Complex<T> centers[9];
            centers[0] = c;
            for (int k = 0; k < 8; ++k) {
                T ang = T(2) * std::acos(T(-1)) * T(k) / T(8);
                centers[k + 1] = c + Complex<T>(T(0.75) * r * std::cos(ang),
                                                T(0.75) * r * std::sin(ang));
            }
            const T rr = r / T(2);
            bool moved = false;
            for (int k = 0; k < 9; ++k) {
                if (hasRootInUnitDisk(shiftScale(work, centers[k], rr), tol)) {
                    c = centers[k]; r = rr; moved = true; break;
                }
            }
            if (!moved) {
                // Fallback: move to the center with the smallest |P| and only shrink the radius
                int best = 0; T bestV = normSq(evalComplex(work, centers[0]));
                for (int k = 1; k < 9; ++k) {
                    T v = normSq(evalComplex(work, centers[k]));
                    if (v < bestV) { bestV = v; best = k; }
                }
                c = centers[best]; r = rr;
            }
        }

        Complex<T> root = newtonPolish(orig, c, 80, tol);   // polish on the original polynomial
        roots.push_back(root);
        work = deflate(work, root);
    }

    // Clean up the imaginary part of real roots
    for (auto& z : roots)
        if (abs(z.im) < tol * (T(1) + abs(z.re))) z.im = T(0);
    return roots;
}

// ================================================================
// B-4. Bairstow method
// ================================================================

/// Polynomial root-finding by the Bairstow method (real-coefficient quadratic factorization)
///
/// Iteratively performs synthetic division of a real-coefficient polynomial by (x² + r·x + s),
/// correcting r, s by Newton's method to extract a quadratic factor.
/// Since the computation stays in real coefficients, complex roots are obtained naturally as conjugate pairs.
///
/// References: 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);

    // Descending-power coefficient array a[0]=highest degree, ..., a[n]=constant
    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) {
        // Initial estimate: based on the ratio of the low-order coefficients
        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) {
            // Synthetic division: divide P(x) by (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];

            // For partial derivatives: apply the same recurrence to b
            // c[j] = b[j] - r·c[j-1] - s·c[j-2]
            // Result: ∂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 system:
            //   [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) {
                // Singular → add a perturbation and retry
                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;
            }
        }

        // Find the roots of the quadratic factor 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)));
        }

        // Deflation: make the quotient polynomial b[0..n-2] the next a
        n -= 2;
        std::vector<T> newA(n + 1);
        for (int i = 0; i <= n; ++i)
            newA[i] = b[i];
        a = std::move(newA);
    }

    // Remaining linear or quadratic
    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. Unified dispatcher
// ================================================================

/// Find all roots of a polynomial (automatically selects the best algorithm by degree)
///
/// Algorithm selection:
///   Degree 1-4: closed-form solution (with cancellation mitigation)
///   Degree 5-19: Jenkins-Traub (one root at a time + deflation, high precision)
///   Degree 20 and above: Aberth-Ehrlich + Newton polish (all roots simultaneously,
///             avoids accumulation of deflation error)
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)
{
    // Remove leading zero coefficients and find the effective degree
    int deg = p.degree();
    if (deg < 0) return {};  // zero polynomial

    // Obtain the effective polynomial (already normalized, so degree() is correct)
    switch (deg) {
        case 0: return {};  // constant polynomial → no roots
        case 1: return solveLinear(p);
        case 2: return solveQuadratic(p);
        case 3: return solveCubic(p);
        case 4: return solveQuartic(p);
        default:
            // Degree 20 and above: DKA (all roots simultaneously + Newton polish)
            // Avoids accumulation of deflation error
            if (deg >= 20)
                return durandKernerAberth(p, eps, maxIter);
            // Degree 5-19: Jenkins-Traub (one root at a time, high precision)
            return jenkinsTraub(p);
    }
}

// ================================================================
// C'. Root-finding with convergence flags ─ *post-verifies* the result to report "did it find all roots / is the precision sufficient"
// ================================================================
//
// Motivation (TODO A.1): solvePolynomial / jenkinsTraub / durandKernerAberth / laguerre / bairstow
//   only return a vector of Complex<T>, and do not tell the caller (a) whether *all* degree-many roots were
//   found (an iteration cutoff may silently return fewer = a counting hole), nor (b) how precise each root is.
//   Without changing the existing signatures at all (non-destructive), add a thin wrapper that *post-verifies* the result.
//   ★ Rather than "trusting the solver's internal converged flag", it actually substitutes the returned roots into the polynomial
//     and measures the residual |p(root)|. Hence even if the solver is buggy and returns too few / sloppy roots, converged will not
//     wrongly become true (the double-side counterpart of numeric_bridge's "don't trust generation, verify it" spirit).
//     If an exact (proven) guarantee is required, also use numeric_bridge::certifyBracketedRoot.
template<typename T>
struct PolyRootsResult {
    std::vector<Complex<T>> roots;       // roots returned by the solver (real roots have im≈0)
    bool   converged   = false;          // all roots found (found==expected) ∧ all residuals ≤ tol
    bool   complete    = false;          // whether degree-many roots were returned (found==expected)
    std::size_t found    = 0;            // number of returned roots
    std::size_t expected = 0;            // expected number of roots (= effective degree)
    T      maxResidual = T(0);           // max_i |p(root_i)|  (over the found roots)
};

/// Find all roots, post-verify the result, and return it with convergence/completeness flags (non-destructive, calls the existing API internally).
///   residualTol: upper bound on |p(root)| (absolute residual). For large-scale polynomials, the caller may
///                judge by its own criterion based on maxResidual (the raw value is also returned).
template<typename T>
PolyRootsResult<T> solvePolynomialChecked(
    const Polynomial<T>& p,
    T residualTol = std::numeric_limits<T>::epsilon() * T(1e4),
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    std::size_t maxIter = 1000)
{
    PolyRootsResult<T> r;
    int deg = p.degree();
    if (deg < 0) return r;                                   // zero polynomial → converged=false (ill-formed)
    r.expected = static_cast<std::size_t>(deg);
    r.roots    = solvePolynomial(p, eps, maxIter);
    r.found    = r.roots.size();
    r.complete = (r.found == r.expected);
    T maxres = T(0);
    for (const auto& z : r.roots) {
        T res = abs(p(z));                                   // |p(root)| (modulus of Complex, ADL)
        if (res > maxres) maxres = res;
    }
    r.maxResidual = maxres;
    r.converged   = r.complete && (deg == 0 || maxres <= residualTol);
    return r;
}

// ================================================================
// D. Real root isolation by Vincent's method (Vincent Real Root Isolation)
// ================================================================

/**
 * @brief Real-root isolation interval
 */
template<typename T>
struct RealRootInterval {
    T lower;           ///< lower bound of the interval
    T upper;           ///< upper bound of the interval
    bool exact;        ///< if true, lower == upper (an exact root)
};

namespace detail {

/// Number of sign variations in a coefficient sequence (Descartes' rule of signs)
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;
}

/// Number of sign variations of a polynomial
template<typename T>
std::size_t signVariations(const Polynomial<T>& p) {
    return signVariations(p.coefficients());
}

/// Taylor shift: compute p(x + c) (Horner scheme, O(n²))
template<typename T>
Polynomial<T> taylorShift(const Polynomial<T>& p, const T& c) {
    int n = p.degree();
    if (n < 0) return p;

    // Transform coeffs[i] into the coefficients of 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));
}

/// Reciprocal polynomial: x^n * p(1/x) (reverse the coefficients)
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 upper bound: upper bound on the positive real roots 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;
}

/// Square-free part: 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;
}

/// Refine a root by bisection
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 the signs are the same, return the midpoint
    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);
}

/// Isolation of positive real roots (VCA bisection)
/// Tracks the Möbius transform 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;  // transformed polynomial
        T a, b, c, d;     // Möbius parameters
    };

    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;

        // Normalize the leading coefficient to be positive
        if (q[deg] < T(0)) q = q * T(-1);

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

        if (v == 1) {
            // Find the positive real root of q and map it back to the original coordinate via the inverse Möbius transform
            T t_root;
            if (deg == 1) {
                // Linear: solve directly
                t_root = -q[0] / q[1];
            } else {
                // Find the positive root of q by bisection
                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;
        }

        // Check for a root at 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});
            // Deflate by 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, ∞) part: 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) part: 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 Real root isolation by Vincent's method
 *
 * Isolates all real roots of a polynomial into mutually disjoint intervals.
 * Uses the bisection version of the Vincent-Collins-Akritas (VCA) algorithm.
 *
 * Each interval contains exactly one real root.
 * When exact == true, it is the exact root value (lower == upper).
 *
 * @param p       input polynomial
 * @param epsilon lower bound tolerance on the interval width (default: 100ε)
 * @return vector of isolation intervals (lower ≤ upper, sorted)
 */
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;

    // Linear: solve directly
    if (deg == 1) {
        T root = -p[0] / p[1];
        result.push_back({root, root, true});
        return result;
    }

    // Extract the root at 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;
    }

    // Square-free part
    auto sf = detail::squareFree(q);
    T bound = detail::cauchyBound(sf);

    // Isolate the positive roots
    detail::vincentPositiveRoots(sf, bound, result, epsilon);

    // Negative roots: isolate the positive roots of p(-x) and flip the sign
    // p(-x): flip the sign of the odd-degree coefficients
    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});
    }

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

    return result;
}

// ================================================================
// Bernoulli method (extraction of the dominant root by the power method)
// ================================================================

/**
 * @brief Find the largest-magnitude root (dominant root) of a polynomial by the Bernoulli method (power method)
 *
 * Make the polynomial p(x) = a_n x^n + ... + a_1 x + a_0 monic, and exploit the property that the ratio
 * y_k / y_{k-1} of the impulse response of the characteristic-equation-based linear recurrence
 * y_k = -c_{n-1} y_{k-1} - ... - c_0 y_{k-n} converges to the dominant root.
 *
 * Reference: Hayato Togawa, "Numerical Computation", p74
 *
 * @param poly polynomial (degree ≥ 1)
 * @param tolerance converges when the difference between successive approximate roots is at most this
 * @param max_iterations maximum number of iterations
 * @return the dominant root and the convergence status
 *
 * @note Does not converge when the largest-magnitude root is a repeated root or an equal-magnitude complex conjugate pair.
 */
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};

    // Make monic: divide by the leading coefficient
    T lead = poly[n];
    std::vector<T> c(n);
    for (int i = 0; i < n; ++i) {
        c[i] = poly[i] / lead;
    }

    // State vector y[0..n-1]: the latest n values of the impulse response
    // Initial values: y[0]=1 (impulse), y[1..n-1]=0
    std::vector<T> y(n, T(0));
    y[0] = T(1);

    // Compute the impulse response via the linear recurrence
    // 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);

    // Skip the first few steps, where y[previous] = 0 can occur
    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];
        }

        // Compute the ratio (after warmup)
        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;
        }

        // Shift: discard y[n-1] and put y_new at the front
        for (int j = n - 1; j > 0; --j) {
            y[j] = y[j - 1];
        }
        y[0] = y_new;

        // Overflow prevention: scale if the values grow too large
        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