Jenkins-Traub from simbody converted into a class template for float/double/... specializations

jenkins-traub-rpoly-qd.cpp

#include <qd/dd_real.h>
#include <qd/qd_real.h>

#include "jenkins-traub-rpoly.ipp"

template struct RPoly<dd_real>;
template struct RPoly<qd_real>;

jenkins-traub-rpoly.cpp

#include <cmath>

using std::max;
using std::abs;
using std::atan;
using std::exp;
using std::log;
using std::pow;
using std::sin;
using std::cos;
using std::sqrt;

namespace {
    int to_int(double arg) { return arg; }
}

#include "jenkins-traub-rpoly.ipp"

template struct RPoly<float>;
template struct RPoly<double>;
template struct RPoly<long double>;

jenkins-traub-rpoly.hpp

#pragma once
#include <memory>
#include <limits>


template <typename T>
struct Complex {
    T re {0}, im {0};
    inline Complex<T> operator+(const T& other) const {
        return Complex<T>{other + re, im};
    }
    inline Complex<T> operator*(const T& other) const {
        return Complex<T>{other*re, other*im};
    }
    inline Complex<T> operator*(const Complex<T>& other) const {
        return Complex<T>{re*other.re - im*other.im, re*other.im + im*other.re};
    }
};

template <typename T>
struct RPoly {
    /*  static constexpr <-- qdlib not constexpr friendly */ T eta {std::numeric_limits<T>::epsilon()};
    const int degree;
private:
    T sr,si,u,v,a,b,c,d,a1,a2;
    T a3,a6,a7,e,f,g,h,szr,szi,lzr,lzi;
    T are,mre;
    int n,nn,nmi,zerok;

    std::unique_ptr<T[]> buffer;
    T * const __restrict temp;
    T * const __restrict pt;
    T * const __restrict p;
    T * const __restrict qp;
    T * const __restrict k;
    T * const __restrict qk;
    T * const __restrict svk;

public:
    RPoly(int degree);
    int findRoots(const T * const coeffs, T * const zeror, T * const zeroi);
    T refineRealRoot(const T * const coeffs, const T& rr) const;
    T evalPoly(const T * const coeffs, T x) const;
    Complex<T> evalPoly(const T * const coeffs, T re, T im) const;
private:
    void quad(T a,T b1,T c,T *sr,T *si,
            T *lr,T *li);
    void fxshfr(int l2, int *nz);
    void quadit(T *uu,T *vv,int *nz);
    void realit(T *sss, int *nz, int *iflag);
    void calcsc(int *type);
    void nextk(int *type);
    void newest(int type,T *uu,T *vv);
    void quadsd(int n,T *u,T *v,T *p,T *q, T *a,T *b);
};

jenkins-traub-rpoly.ipp

#include "jenkins-traub-rpoly.hpp"
/* based on rpoly.cpp from simbody (Apache License 2.0), which is in turns based on TOMS493 (RPOLY) from netlib*/

// #include <summation_cxx/accumulator.hpp>



template <typename T>
RPoly<T>::RPoly(int degree)
    : degree(degree)
    , buffer(std::make_unique<T[]>((degree+1)*7))
    , temp(&buffer[0*(degree+1)])
    , pt(&buffer[1*(degree+1)])
    , p(&buffer[2*(degree+1)])
    , qp(&buffer[3*(degree+1)])
    , k(&buffer[4*(degree+1)])
    , qk(&buffer[5*(degree+1)])
    , svk(&buffer[6*(degree+1)])
{
}

template <typename T>
T RPoly<T>::evalPoly(const T * const coeffs, T x) const {
    T result {0};
    for (int i=0; i<=degree; ++i) { // Horner's method
        result = result*x + coeffs[i];
    }
    return result;
}

template <typename T>
Complex<T> RPoly<T>::evalPoly(const T * const coeffs, T re, T im) const {
    Complex<T> result {};
    Complex<T> x { re, im };
    for (int i=0; i<=degree; ++i) {
        result = result*x + coeffs[i];
    }
    return result;
}

template <typename T>
T RPoly<T>::refineRealRoot(const T *const coeffs, const T &x) const {
    if (degree < 1) {
        return x;
    }
    if (degree > 1) {
        /* single step in Halley's method: dx = -2f(x)*f'(x)/(2f'(x))**2 - f(x)*f''(x)) */
        using Accu_t = T /*summation_cxx::AccumulatorNeumaier<T>*/;
        Accu_t fpp{0}, fp{0}, f{0}, mu{0};
        for (int i=0; i<=degree; ++i) {
            /* Horner's method */
            f = f*x + coeffs[i];
            mu = fabs(x)*mu + fabs(f); /* Algorithm 5.1, p 95, Higham, Accuracy and stability of numerical algo. */
            if (i+1 <= degree) {
                fp = fp*x + coeffs[i]*(degree-i);
            }
            if (i + 2 <= degree) {
                fpp = fpp * x + coeffs[i] * (degree - i) * (degree - i - 1);
            }
        }
        mu = eta*(2*mu - fabs(f)); /* error bound on mu */
        Accu_t dx = -2 * f / (2 * fp - f * fpp / fp);
        T xnext = x + dx/*.template to<T>()*/;
        T fnext = evalPoly(coeffs, xnext);
        if (fnext == 0) {
            return xnext;
        }
        if (fabs(fnext)+mu < fabs(f/*.template to<T>()*/)) {
            return xnext;
        } else {
            return x;
            // xnext = x + 0.5*dx/*.template to<T>()*/;
            // fnext = evalPoly(coeffs, xnext);
            // if (abs(fnext) < abs(f/*.template to<T>()*/)) {
            //     return xnext;
            // } else {
            //     return x;
            // }
        }
    } else {
        /* single step in Newton's method: dx = -f(x)/f'(x) */
        T fp{0}, f{0};
        for (int i=0; i<=degree; ++i) {
            f = f*x + coeffs[i];
            if (i+1 <= degree) {
                fp = fp*x + coeffs[i]*(degree-i);
            }
        }
        T dx = - f / fp;
        return x + dx;
    }

}

template <typename T>
int RPoly<T>::findRoots(const T * const __restrict coeffs_decr, T * const __restrict zeror, T * const __restrict zeroi)
{
    const auto max_repr = std::numeric_limits<T>::max();
    const auto min_repr = std::numeric_limits<T>::min();
    T pi_4 {atan(T(1))};
    T t,aa,bb,cc,factor,rot {T(94)/T(45)*pi_4};
    T lo {min_repr/eta},max,min,xx {sqrt((T) 0.5)};
    T yy {-xx},cosr {cos(rot)},sinr {sin(rot)},xxx,x,sc,bnd;
    T xm,ff,df,dx,base {2.0};
    int cnt,nz,i,j,jj,l,nm1,zerok;

    are = eta;
    mre = eta;
/*  Initialization of constants for shift rotation. */
    n = degree;

/*  Algorithm fails if the leading coefficient is zero (or if degree==0). */
    if (n < 1 || coeffs_decr[0] == 0.0)
        return -1;

/*  Remove the zeros at the origin, if any. */
    while (coeffs_decr[n] == 0.0) {
        j = degree - n;
        zeror[j] = 0.0;
        zeroi[j] = 0.0;
        n--;
    }

    // sherm 20130410: If all coefficients but the leading one were zero, then
    // all solutions are zero; should be a successful (if boring) return.
    if (n == 0)
        return degree;

/*  Make a copy of the coefficients. */
    for (i=0;i<=n;i++)
        p[i] = coeffs_decr[i];
/*  Start the algorithm for one zero. */
_40:
    if (n == 1) {
        zeror[degree-1] = -p[1]/p[0];
        zeroi[degree-1] = 0.0;
        n -= 1;
        goto _99;
    }
/*  Calculate the final zero or pair of zeros. */
    if (n == 2) {
        quad(p[0],p[1],p[2],&zeror[degree-2],&zeroi[degree-2],
            &zeror[degree-1],&zeroi[degree-1]);
        n -= 2;
        goto _99;
    }
/*  Find largest and smallest moduli of coefficients. */
    max = 0.0;
    min = max_repr;
    for (i=0;i<=n;i++) {
        x = fabs(p[i]);
        if (x > max) max = x;
        if (x != 0.0 && x < min) min = x;
    }
/*  Scale if there are large or very small coefficients.
 *  Computes a scale factor to multiply the coefficients of the
 *  polynomial. The scaling si done to avoid overflow and to
 *  avoid undetected underflow interfering with the convergence
 *  criterion. The factor is a power of the base.
 */
    sc = lo/min;
    if (sc > 1.0 && max_repr/sc < max) goto _110;
    if (sc <= 1.0) {
        if (max < 10.0) goto _110;
        if (sc == 0.0)
            sc = min_repr;
    }
    l = to_int((log(sc)/log(base) + 0.5));
    factor = T(pow(T(2.0), l)); // extraneous cast required by VS 16.8.2
    if (factor != 1.0) {
        for (i=0;i<=n;i++)
            p[i] = factor*p[i];     /* Scale polynomial. */
    }
_110:
/*  Compute lower bound on moduli of roots. */
    for (i=0;i<=n;i++) {
        pt[i] = (fabs(p[i]));
    }
    pt[n] = - pt[n];
/*  Compute upper estimate of bound. */
    x = exp((log(-pt[n])-log(pt[0])) / (T)n);
/*  If Newton step at the origin is better, use it. */
    if (pt[n-1] != 0.0) {
        xm = -pt[n]/pt[n-1];
        if (xm < x)  x = xm;
    }
/*  Chop the interval (0,x) until ff <= 0 */
    while (1) {
        xm = x*(T) 0.1;
        ff = pt[0];
        for (i=1;i<=n;i++)
            ff = ff*xm + pt[i];
        if (ff <= 0.0) break;
        x = xm;
    }
    dx = x;
/*  Do Newton interation until x converges to two
 *  decimal places.
 */
    while (fabs(dx/x) > 0.005) {
        ff = pt[0];
        df = ff;
        for (i=1;i<n;i++) {
            ff = ff*x + pt[i];
            df = df*x + ff;
        }
        ff = ff*x + pt[n];
        dx = ff/df;
        x -= dx;
    }
    bnd = x;
/*  Compute the derivative as the initial k polynomial
 *  and do 5 steps with no shift.
 */
    nm1 = n - 1;
    for (i = 1; i < n; i++) {
        k[i] = (T)(n - i) * p[i] / (T)n;
    }
    k[0] = p[0];
    aa = p[n];
    bb = p[n-1];
    zerok = (k[n-1] == 0);
    for(jj=0;jj<5;jj++) {
        cc = k[n-1];
        if (!zerok) {
/*  Use a scaled form of recurrence if value of k at 0 is nonzero. */
            t = -aa/cc;
            for (i=0;i<nm1;i++) {
                j = n-i-1;
                k[j] = t*k[j-1]+p[j];
            }
            k[0] = p[0];
            zerok = (fabs(k[n-1]) <= abs(bb)*eta*10.0);
        }
        else {
/*  Use unscaled form of recurrence. */
            for (i=0;i<nm1;i++) {
                j = n-i-1;
                k[j] = k[j-1];
            }
            k[0] = 0.0;
            zerok = (k[n-1] == 0.0);
        }
    }
/*  Save k for restarts with new shifts. */
    for (i=0;i<n;i++)
        temp[i] = k[i];
/*  Loop to select the quadratic corresponding to each new shift. */
    for (cnt = 0;cnt < 20;cnt++) {
/*  Quadratic corresponds to a double shift to a
 *  non-real point and its complex conjugate. The point
 *  has modulus bnd and amplitude rotated by 94 degrees
 *  from the previous shift.
 */
        xxx = cosr*xx - sinr*yy;
        yy = sinr*xx + cosr*yy;
        xx = xxx;
        sr = bnd*xx;
        si = bnd*yy;
        u = (T) -2.0 * sr;
        v = bnd;
        fxshfr(20*(cnt+1),&nz);
        if (nz != 0) {
/*  The second stage jumps directly to one of the third
 *  stage iterations and returns here if successful.
 *  Deflate the polynomial, store the zero or zeros and
 *  return to the main algorithm.
 */
            j = degree - n;
            zeror[j] = szr;
            zeroi[j] = szi;
            n -= nz;
            for (i=0;i<=n;i++)
                p[i] = qp[i];
            if (nz != 1) {
                zeror[j+1] = lzr;
                zeroi[j+1] = lzi;
            }
            goto _40;
        }
/*  If the iteration is unsuccessful another quadratic
 *  is chosen after restoring k.
 */
        for (i=0;i<n;i++) {
            k[i] = temp[i];
        }
    }
/*  Return with failure if no convergence with 20 shifts. */
_99:
    return degree - n;
}
/*  Computes up to L2 fixed shift k-polynomials,
 *  testing for convergence in the linear or quadratic
 *  case. Initiates one of the variable shift
 *  iterations and returns with the number of zeros
 *  found.
 */
template <typename T>
void RPoly<T>::fxshfr(int l2,int *nz)
{
    T svu,svv,ui,vi,s;
    T betas,betav,oss,ovv,ss,vv,ts,tv;
    T ots,otv,tvv,tss;
    int type, i,j,iflag,vpass,spass,vtry,stry;

    *nz = 0;
    betav = 0.25;
    betas = 0.25;
    oss = sr;
    ovv = v;
/*  Evaluate polynomial by synthetic division. */
    quadsd(n,&u,&v,p,qp,&a,&b);
    calcsc(&type);
    for (j=0;j<l2;j++) {
/*  Calculate next k polynomial and estimate v. */
        nextk(&type);
        calcsc(&type);
        newest(type,&ui,&vi);
        vv = vi;
/*  Estimate s. */
        ss = 0.0;
        if (k[n-1] != 0.0) ss = -p[n]/k[n-1];
        tv = 1.0;
        ts = 1.0;
        if (j == 0 || type == 3) goto _70;
/*  Compute relative measures of convergence of s and v sequences. */
        if (vv != 0.0) tv = fabs((vv-ovv)/vv);
        if (ss != 0.0) ts = fabs((ss-oss)/ss);
/*  If decreasing, multiply two most recent convergence measures. */
        tvv = 1.0;
        if (tv < otv) tvv = tv*otv;
        tss = 1.0;
        if (ts < ots) tss = ts*ots;
/*  Compare with convergence criteria. */
        vpass = (tvv < betav);
        spass = (tss < betas);
        if (!(spass || vpass)) goto _70;
/*  At least one sequence has passed the convergence test.
 *  Store variables before iterating.
 */
        svu = u;
        svv = v;
        for (i=0;i<n;i++) {
            svk[i] = k[i];
        }
        s = ss;
/*  Choose iteration according to the fastest converging
 *  sequence.
 */
        vtry = 0;
        stry = 0;
        if ( ( spass && (!vpass) ) || tss < tvv) goto _40;
_20:
        quadit(&ui,&vi,nz);
        if (*nz > 0) return;
/*  Quadratic iteration has failed. Flag that it has
 *  been tried and decrease the convergence criterion.
 */
        vtry = 1;
        betav *= 0.25;
/*  Try linear iteration if it has not been tried and
 *  the S sequence is converging.
 */
        if (stry || !spass) goto _50;
        for (i=0;i<n;i++) {
            k[i] = svk[i];
        }
_40:
        realit(&s,nz,&iflag);
        if (*nz > 0) return;
/*  Linear iteration has failed. Flag that it has been
 *  tried and decrease the convergence criterion.
 */
        stry = 1;
        betas *=0.25;
        if (iflag == 0) goto _50;
/*  If linear iteration signals an almost double real
 *  zero attempt quadratic iteration.
 */
        ui = -(s+s);
        vi = s*s;
        goto _20;
/*  Restore variables. */
_50:
        u = svu;
        v = svv;
        for (i=0;i<n;i++) {
            k[i] = svk[i];
        }
/*  Try quadratic iteration if it has not been tried
 *  and the V sequence is convergin.
 */
        if (vpass && !vtry) goto _20;
/*  Recompute QP and scalar values to continue the
 *  second stage.
 */
        quadsd(n,&u,&v,p,qp,&a,&b);
        calcsc(&type);
_70:
        ovv = vv;
        oss = ss;
        otv = tv;
        ots = ts;
    }
}
/*  Variable-shift k-polynomial iteration for a
 *  quadratic factor converges only if the zeros are
 *  equimodular or nearly so.
 *  uu, vv - coefficients of starting quadratic.
 *  nz - number of zeros found.
 */
template <typename T>
void RPoly<T>::quadit(T *uu,T *vv,int *nz)
{
    T ui,vi;
    T mp,omp,ee,relstp,t,zm;
    int type,i,j,tried;

    *nz = 0;
    tried = 0;
    u = *uu;
    v = *vv;
    j = 0;
/*  Main loop. */
_10:
    quad(1.0,u,v,&szr,&szi,&lzr,&lzi);
/*  Return if roots of the quadratic are real and not
 *  close to multiple or nearly equal and of opposite
 *  sign.
 */

// sherm 20130410: this early return caused premature termination and
// then failure to find any roots in rare circumstances with 6th order
// ellipsoid nearest point equations. Previously (and in all implementations of
// Jenkins-Traub that I could find) it was just a test on the
// relative size of the difference with respect to the larger root lzr.
// I added the std::max(...,0.1) so that if the roots are small then an
// absolute difference below 0.001 will be considered "close enough" to
// continue iterating. In the case that failed, the roots were around .0001
// and .0002, so their difference was considered large compared with 1% of
// the larger root, 2e-6. But in fact continuing the loop instead resulted in
// six very high-quality roots instead of none.
//
// I'm sorry to say I don't know if I have correctly diagnosed the problem or
// whether this is the right fix! It does pass the regression tests and
// apparently cured the ellipsoid problem. I added the particular bad
// polynomial to the PolynomialTest regression if you want to see it fail.
// These are the problematic coefficients:
//   1.0000000000000000
//   0.021700000000000004
//   2.9889970904696875e-005
//   1.0901272298136685e-008
//  -4.4822782160985054e-012
//  -2.6193432740351220e-015
//  -3.0900602527225053e-019
//
// Original code:
    //if (fabs(fabs(szr)-fabs(lzr)) > 0.01 * fabs(lzr)) return;
// Fixed version:
    if ((T)fabs(fabs(szr) - fabs(lzr)) > (T)0.01 * max((T)fabs(lzr), (T)0.1)) {
        return;
    }
/*  Evaluate polynomial by quadratic synthetic division. */
    quadsd(n,&u,&v,p,qp,&a,&b);
    mp = fabs(a-szr*b) + fabs(szi*b);
/*  Compute a rigorous bound on the rounding error in
 *  evaluating p.
 */
    zm = sqrt(fabs(v));
    ee = (T) 2.0*fabs(qp[0]);
    t = -szr*b;
    for (i=1;i<n;i++) {
        ee = ee*zm + fabs(qp[i]);
    }
    ee = ee*zm + fabs(a+t);
    ee *= ((T)5.0 *mre + (T)4.0*are);
       ee = ee - ((T)5.0*mre+(T)2.0*are)*(fabs(a+t)+fabs(b)*zm)+(T)2.0*are*fabs(t);
/*  Iteration has converged sufficiently if the
 *  polynomial value is less than 20 times this bound.
 */
    if (mp <= 20.0*ee) {
        *nz = 2;
        return;
    }
    j++;
/*  Stop iteration after 20 steps. */
    if (j > 20) return;
    if (j < 2) goto _50;
    if (relstp > 0.01 || mp < omp || tried) goto _50;
/*  A cluster appears to be stalling the convergence.
 *  Five fixed shift steps are taken with a u,v close
 *  to the cluster.
 */
    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 (i=0;i<5;i++) {
        calcsc(&type);
        nextk(&type);
    }
    tried = 1;
    j = 0;
_50:
    omp = mp;
/*  Calculate next k polynomial and new u and v. */
    calcsc(&type);
    nextk(&type);
    calcsc(&type);
    newest(type,&ui,&vi);
/*  If vi is zero the iteration is not converging. */
    if (vi == 0.0) return;
    relstp = fabs((vi-v)/vi);
    u = ui;
    v = vi;
    goto _10;
}
/*  Variable-shift H polynomial iteration for a real zero.
 *  sss - starting iterate
 *  nz  - number of zeros found
 *  iflag - flag to indicate a pair of zeros near real axis.
 */
template <typename T>
void RPoly<T>::realit(T *sss, int *nz, int *iflag)
{
    T pv,kv,t,s;
    T ms,mp,omp,ee;
    int i,j;

    *nz = 0;
    s = *sss;
    *iflag = 0;
    j = 0;
/*  Main loop */
    while (1) {
        pv = p[0];
/*  Evaluate p at s. */
        qp[0] = pv;
        for (i=1;i<=n;i++) {
            pv = pv*s + p[i];
            qp[i] = pv;
        }
        mp = fabs(pv);
/*  Compute a rigorous bound on the error in evaluating p. */
        ms = fabs(s);
        ee = (mre/(are+mre))*fabs(qp[0]);
        for (i=1;i<=n;i++) {
            ee = ee*ms + fabs(qp[i]);
        }
/*  Iteration has converged sufficiently if the polynomial
 *  value is less than 20 times this bound.
 */
        if (mp <= 20.0*((are+mre)*ee-mre*mp)) {
            *nz = 1;
            szr = s;
            szi = 0.0;
            return;
        }
        j++;
/*  Stop iteration after 10 steps. */
        if (j > 10) return;
        if (j < 2) goto _50;
        if (fabs(t) > 0.001*fabs(s-t) || mp < omp) goto _50;
/*  A cluster of zeros near the real axis has been
 *  encountered. Return with iflag set to initiate a
 *  quadratic iteration.
 */
        *iflag = 1;
        *sss = s;
        return;
/*  Return if the polynomial value has increased significantly. */
_50:
        omp = mp;
/*  Compute t, the next polynomial, and the new iterate. */
        kv = k[0];
        qk[0] = kv;
        for (i=1;i<n;i++) {
            kv = kv*s + k[i];
            qk[i] = kv;
        }
        if (fabs(kv) <= fabs(k[n-1])*10.0*eta) {
/*  Use unscaled form. */
            k[0] = 0.0;
            for (i=1;i<n;i++) {
                k[i] = qk[i-1];
            }
        }
        else {
/*  Use the scaled form of the recurrence if the value
 *  of k at s is nonzero.
 */
            t = -pv/kv;
            k[0] = qp[0];
            for (i=1;i<n;i++) {
                k[i] = t*qk[i-1] + qp[i];
            }
        }
        kv = k[0];
        for (i=1;i<n;i++) {
            kv = kv*s + k[i];
        }
        t = 0.0;
        if (fabs(kv) > fabs(k[n-1]*10.0*eta)) t = -pv/kv;
        s += t;
    }
}

/*  This routine calculates scalar quantities used to
 *  compute the next k polynomial and new estimates of
 *  the quadratic coefficients.
 *  type - integer variable set here indicating how the
 *  calculations are normalized to avoid overflow.
 */
template <typename T>
void RPoly<T>::calcsc(int *type)
{
/*  Synthetic division of k by the quadratic 1,u,v */
    quadsd(n-1,&u,&v,k,qk,&c,&d);
    if (fabs(c) > fabs(k[n-1]*100.0*eta)) goto _10;
    if (fabs(d) > fabs(k[n-2]*100.0*eta)) goto _10;
    *type = 3;
/*  Type=3 indicates the quadratic is almost a factor of k. */
    return;
_10:
    if (fabs(d) < fabs(c)) {
        *type = 1;
/*  Type=1 indicates that all formulas are divided by c. */
        e = a/c;
        f = d/c;
        g = u*e;
        h = v*b;
        a3 = a*e + (h/c+g)*b;
        a1 = b - a*(d/c);
        a7 = a + g*d + h*f;
        return;
    }
    *type = 2;
/*  Type=2 indicates that all formulas are divided by d. */
    e = a/d;
    f = c/d;
    g = u*b;
    h = v*b;
    a3 = (a+g)*e + h*(b/d);
    a1 = b*f-a;
    a7 = (f+u)*a + h;
}
/*  Computes the next k polynomials using scalars
 *  computed in calcsc.
 */
template <typename T>
void RPoly<T>::nextk(int *type)
{
    T temp;
    int i;

    if (*type == 3) {
/*  Use unscaled form of the recurrence if type is 3. */
        k[0] = 0.0;
        k[1] = 0.0;
        for (i=2;i<n;i++) {
            k[i] = qk[i-2];
        }
        return;
    }
    temp = a;
    if (*type == 1) {
        temp = b;
    }
    if (fabs(a1) <= fabs(temp)*eta*10.0) {
/*  If a1 is nearly zero then use a special form of the
 *  recurrence.
 */
        k[0] = 0.0;
        k[1] = -a7*qp[0];
        for(i=2;i<n;i++) {
            k[i] = a3*qk[i-2] - a7*qp[i-1];
        }
        return;
    }
/*  Use scaled form of the recurrence. */
    a7 /= a1;
    a3 /= a1;
    k[0] = qp[0];
    k[1] = qp[1] - a7*qp[0];
    for (i=2;i<n;i++) {
        k[i] = a3*qk[i-2] - a7*qp[i-1] + qp[i];
    }
}
/*  Compute new estimates of the quadratic coefficients
 *  using the scalars computed in calcsc.
 */
template <typename T>
void RPoly<T>::newest(int type,T *uu,T *vv)
{
    T a4,a5,b1,b2,c1,c2,c3,c4,temp;

/* Use formulas appropriate to setting of type. */
    if (type == 3) {
/*  If type=3 the quadratic is zeroed. */
        *uu = 0.0;
        *vv = 0.0;
        return;
    }
    if (type == 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;
    }
/*  Evaluate new quadratic coefficients. */
    b1 = -k[n-1]/p[n];
    b2 = -(k[n-2]+b1*p[n-1])/p[n];
    c1 = v*b2*a1;
    c2 = b1*a7;
    c3 = b1*b1*a3;
    c4 = c1 - c2 - c3;
    temp = a5 + b1*a4 - c4;
    if (temp == 0.0) {
        *uu = 0.0;
        *vv = 0.0;
        return;
    }
    *uu = u - (u*(c3+c2)+v*(b1*a1+b2*a7))/temp;
    *vv = v*((T)1.0+c4/temp);
    return;
}

/*  Divides p by the quadratic 1,u,v placing the quotient
 *  in q and the remainder in a,b.
 */
template <typename T>
void RPoly<T>::quadsd(int nn,T *u,T *v,T *p,T *q,
    T *a,T *b)
{
    T c;
    int i;
    *b = p[0];
    q[0] = *b;
    *a = p[1] - (*b)*(*u);
    q[1] = *a;
    for (i=2;i<=nn;i++) {
        c = p[i] - (*a)*(*u) - (*b)*(*v);
        q[i] = c;
        *b = *a;
        *a = c;
    }
}
/*  Calculate the zeros of the quadratic a*z^2 + b1*z + c.
 *  The quadratic formula, modified to avoid overflow, is used
 *  to find the larger zero if the zeros are real and both
 *  are complex. The smaller real zero is found directly from
 *  the product of the zeros c/a.
 */
template <typename T>
void RPoly<T>::quad(T a,T b1,T c,T *sr,T *si,
        T *lr,T *li)
{
        T b,d,e;

        if (a == 0.0) {         /* less than two roots */
            if (b1 != 0.0) {
                *sr = -c / b1;
            } else {
                *sr = 0.0;
            }
            *lr = 0.0;
            *si = 0.0;
            *li = 0.0;
            return;
        }
        if (c == 0.0) {         /* one real root, one zero root */
            *sr = 0.0;
            *lr = -b1/a;
            *si = 0.0;
            *li = 0.0;
            return;
        }
/* Compute discriminant avoiding overflow. */
        b = b1/(T) 2.0;
        if (fabs(b) < fabs(c)) {
            if (c < 0.0) {
                e = -a;
            } else {
                e = a;
            }
            e = b*(b/fabs(c)) - e;
            d = sqrt(fabs(e))*sqrt(fabs(c));
        }
        else {
            e = (T) 1.0 - (a/b)*(c/b);
            d = sqrt(fabs(e))*fabs(b);
        }
        if (e < 0.0) {      /* complex conjugate zeros */
            *sr = -b/a;
            *lr = *sr;
            *si = fabs(d/a);
            *li = -(*si);
        }
        else {
            if (b >= 0.0) { /* real zeros. */
                d = -d;
            }
            *lr = (-b+d)/a;
            *sr = 0.0;
            if (*lr != 0.0) {
                *sr = (c / *lr) / a;
            }
            *si = 0.0;
            *li = 0.0;
        }
}

LICENSE.txt

Apache License
Version 2.0, January 2004
http://www.apache.org/licenses/

TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION

1. Definitions.

"License" shall mean the terms and conditions for use, reproduction, and
distribution as defined by Sections 1 through 9 of this document.

"Licensor" shall mean the copyright owner or entity authorized by the copyright
owner that is granting the License.

"Legal Entity" shall mean the union of the acting entity and all other entities
that control, are controlled by, or are under common control with that entity.
For the purposes of this definition, "control" means (i) the power, direct or
indirect, to cause the direction or management of such entity, whether by
contract or otherwise, or (ii) ownership of fifty percent (50%) or more of the
outstanding shares, or (iii) beneficial ownership of such entity.

"You" (or "Your") shall mean an individual or Legal Entity exercising
permissions granted by this License.

"Source" form shall mean the preferred form for making modifications, including
but not limited to software source code, documentation source, and configuration
files.

"Object" form shall mean any form resulting from mechanical transformation or
translation of a Source form, including but not limited to compiled object code,
generated documentation, and conversions to other media types.

"Work" shall mean the work of authorship, whether in Source or Object form, made
available under the License, as indicated by a copyright notice that is included
in or attached to the work (an example is provided in the Appendix below).

"Derivative Works" shall mean any work, whether in Source or Object form, that
is based on (or derived from) the Work and for which the editorial revisions,
annotations, elaborations, or other modifications represent, as a whole, an
original work of authorship. For the purposes of this License, Derivative Works
shall not include works that remain separable from, or merely link (or bind by
name) to the interfaces of, the Work and Derivative Works thereof.

"Contribution" shall mean any work of authorship, including the original version
of the Work and any modifications or additions to that Work or Derivative Works
thereof, that is intentionally submitted to Licensor for inclusion in the Work
by the copyright owner or by an individual or Legal Entity authorized to submit
on behalf of the copyright owner. For the purposes of this definition,
"submitted" means any form of electronic, verbal, or written communication sent
to the Licensor or its representatives, including but not limited to
communication on electronic mailing lists, source code control systems, and
issue tracking systems that are managed by, or on behalf of, the Licensor for
the purpose of discussing and improving the Work, but excluding communication
that is conspicuously marked or otherwise designated in writing by the copyright
owner as "Not a Contribution."

"Contributor" shall mean Licensor and any individual or Legal Entity on behalf
of whom a Contribution has been received by Licensor and subsequently
incorporated within the Work.

2. Grant of Copyright License.

Subject to the terms and conditions of this License, each Contributor hereby
grants to You a perpetual, worldwide, non-exclusive, no-charge, royalty-free,
irrevocable copyright license to reproduce, prepare Derivative Works of,
publicly display, publicly perform, sublicense, and distribute the Work and such
Derivative Works in Source or Object form.

3. Grant of Patent License.

Subject to the terms and conditions of this License, each Contributor hereby
grants to You a perpetual, worldwide, non-exclusive, no-charge, royalty-free,
irrevocable (except as stated in this section) patent license to make, have
made, use, offer to sell, sell, import, and otherwise transfer the Work, where
such license applies only to those patent claims licensable by such Contributor
that are necessarily infringed by their Contribution(s) alone or by combination
of their Contribution(s) with the Work to which such Contribution(s) was
submitted. If You institute patent litigation against any entity (including a
cross-claim or counterclaim in a lawsuit) alleging that the Work or a
Contribution incorporated within the Work constitutes direct or contributory
patent infringement, then any patent licenses granted to You under this License
for that Work shall terminate as of the date such litigation is filed.

4. Redistribution.

You may reproduce and distribute copies of the Work or Derivative Works thereof
in any medium, with or without modifications, and in Source or Object form,
provided that You meet the following conditions:

You must give any other recipients of the Work or Derivative Works a copy of
this License; and
You must cause any modified files to carry prominent notices stating that You
changed the files; and
You must retain, in the Source form of any Derivative Works that You distribute,
all copyright, patent, trademark, and attribution notices from the Source form
of the Work, excluding those notices that do not pertain to any part of the
Derivative Works; and
If the Work includes a "NOTICE" text file as part of its distribution, then any
Derivative Works that You distribute must include a readable copy of the
attribution notices contained within such NOTICE file, excluding those notices
that do not pertain to any part of the Derivative Works, in at least one of the
following places: within a NOTICE text file distributed as part of the
Derivative Works; within the Source form or documentation, if provided along
with the Derivative Works; or, within a display generated by the Derivative
Works, if and wherever such third-party notices normally appear. The contents of
the NOTICE file are for informational purposes only and do not modify the
License. You may add Your own attribution notices within Derivative Works that
You distribute, alongside or as an addendum to the NOTICE text from the Work,
provided that such additional attribution notices cannot be construed as
modifying the License.
You may add Your own copyright statement to Your modifications and may provide
additional or different license terms and conditions for use, reproduction, or
distribution of Your modifications, or for any such Derivative Works as a whole,
provided Your use, reproduction, and distribution of the Work otherwise complies
with the conditions stated in this License.

5. Submission of Contributions.

Unless You explicitly state otherwise, any Contribution intentionally submitted
for inclusion in the Work by You to the Licensor shall be under the terms and
conditions of this License, without any additional terms or conditions.
Notwithstanding the above, nothing herein shall supersede or modify the terms of
any separate license agreement you may have executed with Licensor regarding
such Contributions.

6. Trademarks.

This License does not grant permission to use the trade names, trademarks,
service marks, or product names of the Licensor, except as required for
reasonable and customary use in describing the origin of the Work and
reproducing the content of the NOTICE file.

7. Disclaimer of Warranty.

Unless required by applicable law or agreed to in writing, Licensor provides the
Work (and each Contributor provides its Contributions) on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied,
including, without limitation, any warranties or conditions of TITLE,
NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A PARTICULAR PURPOSE. You are
solely responsible for determining the appropriateness of using or
redistributing the Work and assume any risks associated with Your exercise of
permissions under this License.

8. Limitation of Liability.

In no event and under no legal theory, whether in tort (including negligence),
contract, or otherwise, unless required by applicable law (such as deliberate
and grossly negligent acts) or agreed to in writing, shall any Contributor be
liable to You for damages, including any direct, indirect, special, incidental,
or consequential damages of any character arising as a result of this License or
out of the use or inability to use the Work (including but not limited to
damages for loss of goodwill, work stoppage, computer failure or malfunction, or
any and all other commercial damages or losses), even if such Contributor has
been advised of the possibility of such damages.

9. Accepting Warranty or Additional Liability.

While redistributing the Work or Derivative Works thereof, You may choose to
offer, and charge a fee for, acceptance of support, warranty, indemnity, or
other liability obligations and/or rights consistent with this License. However,
in accepting such obligations, You may act only on Your own behalf and on Your
sole responsibility, not on behalf of any other Contributor, and only if You
agree to indemnify, defend, and hold each Contributor harmless for any liability
incurred by, or claims asserted against, such Contributor by reason of your
accepting any such warranty or additional liability.

END OF TERMS AND CONDITIONS

APPENDIX: How to apply the Apache License to your work

To apply the Apache License to your work, attach the following boilerplate
notice, with the fields enclosed by brackets "[]" replaced with your own
identifying information. (Don't include the brackets!) The text should be
enclosed in the appropriate comment syntax for the file format. We also
recommend that a file or class name and description of purpose be included on
the same "printed page" as the copyright notice for easier identification within
third-party archives.

   Copyright [yyyy] [name of copyright owner]

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.

Makefile

SHELL = bash
MAKEFILE_DIR:=$(strip $(shell dirname $(realpath $(lastword $(MAKEFILE_LIST)))))
FLAGS_WARN ?=\
 -Wall \
 -Wextra \
 -pedantic \
 -Wdouble-promotion \
 -Wformat=2 \
 -Winvalid-pch \
 -Wnull-dereference \
 -Wmissing-include-dirs \
 -Wswitch-enum \
 -Wswitch-default \
 -Wshadow
ifeq (,$(findstring clang++,$(CXX)))
	FLAGS_WARN:=$(FLAGS_WARN) \
 -Wduplicated-branches \
 -Wduplicated-cond \
 -Wlogical-op \
 -Wrestrict \
 -Wuseless-cast
endif

FAILURE_HANDLER=exit 1; echo 

SRC ?=$(MAKEFILE_DIR)
OPTFLAG ?=-O0
DEBNFO ?=-g3
CXXLD ?= $(CXX)
AR ?= ar # llvm-ar
CONTEXT ?= # e.g. valgrind, /usr/bin/time -v, LSAN_OPTIONS=detect_leaks=0 gdb -q -ex r -args

DEFINES := 
CPPFLAGS :=-I$(MAKEFILE_DIR)/.. -I$(SRC)
CXXFLAGS :=$(FLAGS_WARN) -Werror -std=c++20 $(OPTFLAG) $(DEBNFO) -fPIC -fvisibility=hidden $(CXXFLAGS) $(DEFINES)
LDFLAGS := $(shell eval echo $$LDFLAGS)

# reference: http://make.mad-scientist.net/papers/advanced-auto-dependency-generation/#tldr
TEST_EXES := \
 test-jenkins-traub-rpoly \
 test-jenkins-traub-rpoly-qd \
 test-jenkins-traub-rpoly-boostmp


.PHONY: test clean

TEST_DEPS := $(TEST_EXES)

test: $(TEST_DEPS)
	$(CONTEXT) ./test-jenkins-traub-rpoly
	$(CONTEXT) ./test-jenkins-traub-rpoly-qd
	$(CONTEXT) ./test-jenkins-traub-rpoly-boostmp

clean:
	$(RM) $(TEST_DEPS)
	$(RM) -r ./build/temp.*/ ./build/lib.*/

test-jenkins-traub-rpoly: test-jenkins-traub-rpoly.o jenkins-traub-rpoly.o 
	$(CXXLD) -o $@ $(CXXFLAGS) $^ $(LDFLAGS)

jenkins-traub-rpoly-qd.o: $(SRC)/jenkins-traub-rpoly-qd.cpp jenkins-traub-rpoly.ipp jenkins-traub-rpoly.hpp
	$(CXX) -c -o $@ $(CXXFLAGS) $(QDLIB_CXXFLAGS) $< $(LDFLAGS)

test-jenkins-traub-rpoly-qd: $(SRC)/test-jenkins-traub-rpoly-qd.cpp jenkins-traub-rpoly-qd.o 
	$(CXXLD) -o $@ $(CXXFLAGS) $(QDLIB_CXXFLAGS) $^ $(LDFLAGS) $(LDLIBS) $(QDLIB_LDFLAGS) -lqd

test-jenkins-traub-rpoly-boostmp: test-jenkins-traub-rpoly-boostmp.o
	$(CXXLD) -o $@ $(CXXFLAGS) $^ $(LDFLAGS)

test-jenkins-traub-rpoly-boostmp.cpp

#include <boost/multiprecision/cpp_dec_float.hpp>
namespace {
    int to_int(boost::multiprecision::cpp_dec_float_50 arg) { return static_cast<int>(arg); }
}
// #ifndef __clang__
// // https://github.com/boostorg/math/issues/181#issuecomment-597061742
// #include <boost/multiprecision/float128.hpp>
// namespace {
//     int to_int(boost::multiprecision::float128 arg) { return static_cast<int>(arg); }
// }
// #endif
#include "jenkins-traub-rpoly.hpp"
#include "jenkins-traub-rpoly.ipp"
#include <algorithm>
#include "testing.hpp"

#include "test-jenkins-traub-rpoly.ipp"


int main() {
    int result = 0;
    result += test_rpoly<boost::multiprecision::cpp_dec_float_50>(1.0);
    result += test_rpoly_from_paper<boost::multiprecision::cpp_dec_float_50>(4e38 /* TODO: outrageous */);
    result += test_rpoly_ill_conditioned3a<boost::multiprecision::cpp_dec_float_50>(2e47 /* TODO: outrageous */);
// #ifndef __clang__
//     result += test_rpoly<boost::multiprecision::float128>();
// #endif
    return result;
}

test-jenkins-traub-rpoly-qd.cpp

#include "jenkins-traub-rpoly.hpp"
#include <qd/dd_real.h>
#include <qd/qd_real.h>
#include <algorithm>
#include "testing.hpp"

#include "test-jenkins-traub-rpoly.ipp"

int main() {
    int result = 0;
    result += test_rpoly<dd_real>(6.0);
    result += test_rpoly<qd_real>(6.0);
    result += test_rpoly_from_paper<dd_real>(8e20 /* TODO: outrageous */);
    result += test_rpoly_from_paper<qd_real>(4e52 /* TODO: outrageous */);
    result += test_rpoly_ill_conditioned3a<dd_real>(1.0);
    result += test_rpoly_ill_conditioned3a<qd_real>(1.0);
    return result;
}

test-jenkins-traub-rpoly.cpp

#include "jenkins-traub-rpoly.hpp"
#include "testing.hpp"

using std::abs;
using std::max;
#include <cmath>
using std::log;
#include "test-jenkins-traub-rpoly.ipp"


int main() {
    int result = 0;
    result += test_rpoly<double>(1.0);
    result += test_rpoly<float>(30.0);
    result += test_rpoly<long double>(5.0);
    result += test_rpoly_from_paper<double>(1.0);
    result += test_rpoly_from_paper<float>(30.0);
    result += test_rpoly_from_paper<long double>(4e8 /* TODO: outrageous, fix RPoly. */);
    result += test_rpoly_ill_conditioned3a<double>(1.0);
    result += test_rpoly_ill_conditioned3a<float>(1.0);
    result += test_rpoly_ill_conditioned3a<long double>(1.0);
    return result;
}

test-jenkins-traub-rpoly.ipp

#include "jenkins-traub-rpoly.hpp"  /* for IDE */
#include "testing.hpp" /* for IDE */


// https://github.com/tab58/jenkins-traub-rpoly/blob/master/test.js

// https://github.com/simbody/simbody/blob/a8f49c84e98ccf3b7e6f05db55a29520e5f9c176/SimTKcommon/tests/PolynomialTest.cpp

namespace {
    template <typename T>
    void ignore(T&) { }

    template <typename T>
    int verify_roots(const RPoly<T>& rp, const T * const coeffs, const T * const roots_r, const T * const roots_i, const T atol) {
        for (int i=0; i<rp.degree; ++i) {
            const Complex<T> ev = rp.evalPoly(coeffs, roots_r[i], roots_i[i]);
            if (fabs(ev.re) > atol || fabs(ev.im) > atol) {
                return i;
            }
        }
        return -1;
    }
    template <typename T>
    int verify_real_roots(const RPoly<T>& rp, const T * const coeffs, const T * const roots_r, const T * const roots_i, const T atol) {
        for (int i=0; i<rp.degree; ++i) {
            if (fabs(roots_i[i]) > atol) {
                return i;
            }
            const T ev = rp.evalPoly(coeffs, roots_r[i]);
            if (fabs(ev) > atol) {
                return i;
            }
        }
        return -1;
    }

    template <typename T>
    int check_roots(const int degree, const T * const roots_r, const T * const roots_i, const T * const ref_r, const T * const ref_i, const T atol, const T rtol) {
        assert(rtol >= 0);
        assert(atol >= 0);
        for (int i=0; i<degree; ++i) {
            T delta_r = roots_r[i] - ref_r[i];
            T delta_i = roots_i[i] - ref_i[i];
            double ddelta = delta_r.template convert_to<double>();
            double overshoot = (fabs(delta_r) / (fabs(ref_r[i]) * rtol + atol)).template convert_to<double>();
            ignore(overshoot); ignore(ddelta);
            if (fabs(ref_r[i]) * rtol + atol < fabs(delta_r)) {
                return i;
            }
            if (fabs(ref_i[i]) * rtol + atol < fabs(delta_i)) {
                return i;
            }
        }
        return -1;
    }
}

template <typename T>
int test_rpoly_ill_conditioned3a(double slack) {
    const T eps = std::numeric_limits<T>::epsilon();
    const T N = exp2(T(std::numeric_limits<T>::digits - 7));
    T coeffs[4] = {N+1, 1 - N, -1 - N, N - 1};
    RPoly<T> rp3(3);
    T re3[3], im3[3];
    int nroots = rp3.findRoots(coeffs, re3, im3);
    if (nroots != 3) {
        RETURN_(1);
    }
    if (verify_roots<T>(rp3, coeffs, re3, im3, 2.1) >= 0) {
        RETURN_(1);
    }
    const T re_ref[3] = {1, -1, (N-1)/(N+1)};
    const T im_ref[3] = {0, 0, 0};
    if (check_roots<T>(3, re3, im3, re_ref, im_ref, 100*eps*slack, 100*eps*slack) >= 0){
        RETURN_(1);
    }
    RETURN_(0);
}

template <typename T>
int test_rpoly_from_paper(double slack)
{
    const T eps = std::numeric_limits<T>::epsilon();
    RPoly<T> rp7(7);
    T re7[7], im7[7];
    {
        T coeffs[8] = {1, -6.01, 12.54, -8.545, -5.505, 12.545, -8.035, 2.01};
        int nroots = rp7.findRoots(coeffs, re7, im7);
        if (nroots != 7) {
            RETURN_(1);
        }
        /*
          In [1]: from flint import *
          In [4]: p = fmpq_poly([2010, -8035, 12545, -5505, -8545, 12540, -6010, 1000], 1000)
          In [5]: p.roots()
          Out[5]:
          [([-1.00000000000000 +/- 1e-19], 1),
          ([2.00000000000000 +/- 1e-19], 1),
          ([2.01000000000000 +/- 1e-19], 1),
          (0.500000000000000 + 0.500000000000000j, 1),
          (0.500000000000000 - 0.500000000000000j, 1),
          (1.00000000000000, 2)]
         */
        const T re_ref[7] = { 0.5, 0.5, -1.0, 1.0, 1.0, 2.0, 2.01 };
        const T im_ref[7] = { 0.5, -0.5, 0, 0, 0, 0, 0};
        T re7r[7];
        for (int i=0; i<7; ++i) {
            T err = re7[i] - re_ref[i];
            if (im_ref[i] == 0) {
                T refined = rp7.refineRealRoot(coeffs, re7[i]);
                T err_refined  = refined - re_ref[i];
                if (fabs(err_refined) > (1 + slack/67)*fabs(err)) {
                    RETURN_(1);
                }
                re7r[i] = refined;
            } else {
                re7r[i] = re7[i];
            }
        }
        if (verify_roots<T>(rp7, coeffs, re7r, im7, 3e4*eps*slack) >= 0) {
            RETURN_(1);
        }
        T atol = T(1000*eps) > T(1e-9) ? T(1000*eps) : T(1e-9) /* 1e-9 here should be based on eps... */;
        if (check_roots<T>(7, re7r, im7, re_ref, im_ref, atol, 1000*eps*slack) >= 0){
            RETURN_(1);
        }
    }
    RETURN_(0);
}
template <typename T>
int test_rpoly(double slack) {
    const T eps = std::numeric_limits<T>::epsilon();
    RPoly<T> rp1(1);
    {
        T coeffs[2] = {2, 3}; /* 2*x + 3 */
        T re, im;
        int nroots = rp1.findRoots(coeffs, &re, &im);
        if (nroots != 1) {
            RETURN_(1);
        }
        const T half { 0.5 };
        if (fabs(re + 3*half) > 3*eps) {
            RETURN_(1);
        }
    }
    RPoly<T> rp2(2);
    {
        T work[3+2+2] = {-3, 2, 5, -99, -99, -99, -99}; /* -3x**2 + 2*x + 5*/
        int nroots = rp2.findRoots(work, work + 3, work + 5);
        if (nroots != 2) {
            RETURN_(2);
        }
        const T ref_roots_r[2] = {-1, T(5)/3};
        const T ref_roots_i[2] = {0, 0};
        if (verify_real_roots<T>(rp2, work, work+3, work+5, 10*eps*slack) != -1) {
            RETURN_(1);
        }
        if (check_roots<T>(2, work+3, work+5, ref_roots_r, ref_roots_i, 5*eps*slack, 5*eps*slack) != -1) {
            RETURN_(1);
        }
    }
    RPoly<T> rp4(4);
    T zeror4[4];
    T zeroi4[4];
    {
        // >>> ((x-4)*(x-3)*(x-2)*(x-1)).expand()
        //    4       3       2
        //   x  - 10⋅x  + 35⋅x  - 50⋅x + 24
        //
        T coeffs[5] = {1, -10, 35, -50, 24};
        int nroots = rp4.findRoots(coeffs, zeror4, zeroi4);
        if (nroots != 4) {
            RETURN_(1);
        }
        const T ref_roots_r[4] = {1, 2, 3, 4};
        const T ref_roots_i[4] = {0, 0, 0, 0};
        if (verify_real_roots<T>(rp4, coeffs, zeror4, zeroi4, 250*eps*slack) >= 0) {
            RETURN_(1);
        }
        if (check_roots<T>(rp4.degree, zeror4, zeroi4, ref_roots_r, ref_roots_i, 10*eps*slack, 10*eps*slack) >= 0) {
            RETURN_(1);
        }
    }
    {
        // In [1]: 24*x**4 - 50*x**3 + 35*x**2 -10*x + 1
        // Out[1]:
        //     4       3       2
        // 24⋅x  - 50⋅x  + 35⋅x  - 10⋅x + 1

        // In [2]: solveset(24*x**4 - 50*x**3 + 35*x**2 -10*x + 1)
        // Out[2]: {1/4, 1/3, 1/2, 1}
        T coeffs[5] = {24, -50, 35, -10, 1};
        int nroots = rp4.findRoots(coeffs, zeror4, zeroi4);
        if (nroots != 4) {
            RETURN_(1);
        }
        const T ref_roots_r[4] = { 1/T(4), 1/T(3), 1/T(2), 1 };
        const T ref_roots_i[4] = {0, 0, 0, 0};
        if (check_roots<T>(4, zeror4, zeroi4, ref_roots_r, ref_roots_i, 10*eps*slack, 10*eps*slack) >= 0) {
            RETURN_(1);
        }
    }
    RPoly<T> rp6(6);
    T re6[6], im6[6];
    {
        T coeffs[7] = { 1., 0.093099999999999988, 0.00057211940879762883, 1.4343090324181468e-6,
            1.6097307763625053e-9, 6.6189348690786845e-13, -3.0418139616145048e-18 };

        int nroots = rp6.findRoots(coeffs, re6, im6);
        if (nroots != 6) {
            RETURN_(1);
        }
        if (verify_roots<T>(rp6, coeffs, re6, im6, 10*eps*slack) >= 0) {
            RETURN_(1);
        }
        const T re_ref[6] = {4.54517869625002e-6, -0.00122499999999999, -0.00122499999999999 , -0.00198292148034529,
            -0.00198292148034529, -0.0866887022180057 };
        const T im_ref[6] = {0, 1.12651985673021e-10, -1.12651985673021e-10, 0.0011011665837654, -0.0011011665837654, 0};
        if (check_roots<T>(6, re6, im6, re_ref, im_ref, 1.2e-10, 10*eps*slack) >= 0){
            RETURN_(1);
        }
    }
    {
        T coeffs[7] = { 1.,  0.021700000000000004,  0.00013355532616269355,  1.8068473626980300e-7,  6.2414644021223684e-11, 6.4036661086410838e-15, -6.8627441483352105e-22 };
        int nroots = rp6.findRoots(coeffs, re6, im6);
        if (nroots != 6) {
            RETURN_(1);
        }
        if (verify_roots<T>(rp6, coeffs, re6, im6, 10*eps*slack) >= 0) {
            RETURN_(1);
        }
        const T re_ref[6] = { 1.07057243683631e-7, -0.000226941898625434 , -0.000226941898625434,
            -0.00124724614073959, -0.00896457624670489, -0.0110344008725483};
        const T im_ref[6] = { 0.0, +2.13356769572618e-5, -2.13356769572618e-5, 0.0, 0.0, 0.0};
        if (check_roots<T>(6, re6, im6, re_ref, im_ref, 1.2e-10, 10*eps*slack) >= 0){
            RETURN_(1);
        }
    }
    {
        T coeffs[7] = {  1.,  0.021700000000000004,  2.9889970904696875e-5,  1.0901272298136685e-8, -4.4822782160985054e-12, -2.6193432740351220e-15, -3.0900602527225053e-19};
        int nroots = rp6.findRoots(coeffs, re6, im6);
        if (nroots != 6) {
            RETURN_(1);
        }
        if (verify_roots<T>(rp6, coeffs, re6, im6, 10*eps*slack) >= 0) {
            RETURN_(1);
        }
        const T re_ref[6] = { -0.000227633733894355 , -0.000227633733894355, 0.000433826366836664, -0.000713709028258555, -0.000713709028258555, -0.0202511408425308};
        const T im_ref[6] = { +3.53990853672026e-5, -3.53990853672026e-5, 0.0, +0.000391625935770687, -0.000391625935770687, 0.0 };
        if (check_roots<T>(6, re6, im6, re_ref, im_ref, 1.2e-10, 10*eps*slack) >= 0){
            RETURN_(1);
        }
    }
    RETURN_(0);
}