PolynomialRoots.hh

Go to the documentation of this file.

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2014                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Università degli Studi di Trento                                    |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/
 
#ifndef POLYNOMIAL_ROOTS_HH
#define POLYNOMIAL_ROOTS_HH
 
#include <cmath>
#include <cfloat>
#include <iostream>
 
#include <complex>
#include <cstdint>
#include <format>
#include <iomanip>
#include <iosfwd>
#include <limits>
#include <sstream>
#include <string>

 
#if POLYNOMIAL_ROOTS_HAS_MULTIPRECISION
  #include <boost/multiprecision/cpp_bin_float.hpp>
  #include <boost/math/special_functions/cbrt.hpp>
#endif
 
namespace PolynomialRoots
{

  using real_type    = double;
  using integer      = int;
  using ostream_type = std::basic_ostream<char>;
  using istream_type = std::basic_istream<char>;

  inline constexpr integer MAXDEGREE = 100;

  template <typename... Args> inline void root_assert( bool cond, std::format_string<Args...> fmt, Args &&... args )
  {
    if ( !cond ) std::runtime_error( std::format( fmt, std::forward<Args>( args )... ) );
  }
 
#if POLYNOMIAL_ROOTS_HAS_MULTIPRECISION
  using quad_real = boost::multiprecision::cpp_bin_float_100;
#endif
 
}  // namespace PolynomialRoots
 
namespace std
{
 
#if POLYNOMIAL_ROOTS_HAS_MULTIPRECISION
  template <>
 
  struct formatter<PolynomialRoots::quad_real, char> : formatter<std::string, char>
 
  {
    auto format( PolynomialRoots::quad_real const & x, format_context & ctx ) const
    {
      std::ostringstream ss;
      constexpr int      digits = std::numeric_limits<PolynomialRoots::quad_real>::max_digits10 > 0
                                    ? std::numeric_limits<PolynomialRoots::quad_real>::max_digits10
                                    : 36;
      ss << std::setprecision( digits ) << x;
      return formatter<std::string, char>::format( ss.str(), ctx );
    }
  };
#endif
}  // namespace std
 
#include "PolynomialRoots-complex.hxx"
 
namespace PolynomialRoots
{
  using std::isfinite;

  template <typename T_real> T_real machepsiT();
  template <typename T_real> T_real toleranceT();

  template <typename T_real> T_real eval_poly( T_real const op[], integer Degree, T_real const & x );

  template <typename T_real>
  void eval_poly_Dpoly( T_real const op[], integer Degree, T_real const & x, T_real & p, T_real & dp );

  template <typename T_real> bool Newton_step( T_real const op[], integer Degree, T_real & x );

  template <typename T_real, typename T_complex>
  T_complex eval_poly_complex( T_real const op[], integer Degree, T_complex const & x );

  [[nodiscard]] int roots( real_type const * op, integer Degree, real_type * zeror, real_type * zeroi );
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  /*\
   |    ___                  _           _   _
   |   / _ \ _   _  __ _  __| |_ __ __ _| |_(_) ___
   |  | | | | | | |/ _` |/ _` | '__/ _` | __| |/ __|
   |  | |_| | |_| | (_| | (_| | | | (_| | |_| | (__
   |   \__\_\\__,_|\__,_|\__,_|_|  \__,_|\__|_|\___|
   |
   |  A * x^2 + B * x + C
  \*/
  template <typename T_real, typename T_complex> class QuadraticT
  {
    T_real  m_ABC[3]{ 0, 0, 0 };
    T_real  m_r0   = 0;
    T_real  m_r1   = 0;
    integer m_nrts = 0;
    bool    m_cplx = false;
    bool    m_dblx = false;
 
    void find_roots();
 
  public:
    using value_type   = T_real;
    using complex_type = T_complex;

    QuadraticT() = default;

    QuadraticT( T_real const & a, T_real const & b, T_real const & c )
    {
      m_ABC[0] = a;
      m_ABC[1] = b;
      m_ABC[2] = c;
      // find roots only on finite values
      root_assert(
        isfinite( a ) && isfinite( b ) && isfinite( c ),
        "QuadraticT( a={}, b={}, c={} ) arguments must be finite!",
        a, b, c
      );
      find_roots();
    }

    void setup( T_real const & a, T_real const & b, T_real const & c )
    {
      m_ABC[0] = a;
      m_ABC[1] = b;
      m_ABC[2] = c;
      find_roots();
    }

    integer num_roots() const { return m_nrts; }

    integer numRoots() const { return m_nrts; }

    bool complex_root() const { return m_cplx; }

    bool complexRoot() const { return m_cplx; }

    bool double_root() const { return m_dblx; }

    bool doubleRoot() const { return m_dblx; }

    integer get_real_roots( T_real r[] ) const;

    integer getRealRoots( T_real r[] ) const { return get_real_roots( r ); }

    integer get_positive_roots( T_real r[] ) const;

    integer getPositiveRoots( T_real r[] ) const { return get_positive_roots( r ); }

    integer get_negative_roots( T_real r[] ) const;

    integer getNegativeRoots( T_real r[] ) const { return get_negative_roots( r ); }

    integer get_roots_in_range( T_real const & a, T_real const & b, T_real r[] ) const;

    integer getRootsInRange( T_real const & a, T_real const & b, T_real r[] ) const
    { return get_roots_in_range( a, b, r ); }

    integer get_roots_in_open_range( T_real const & a, T_real const & b, T_real r[] ) const;

    integer getRootsInOpenRange( T_real const & a, T_real const & b, T_real r[] ) const
    { return get_roots_in_open_range( a, b, r ); }

    T_real real_root0() const { return m_r0; }

    T_real real_root1() const { return m_r1; }

    T_complex root0() const { return m_cplx ? T_complex( m_r0, m_r1 ) : T_complex( m_r0, 0 ); }

    T_complex root1() const { return m_cplx ? T_complex( m_r0, -m_r1 ) : T_complex( m_r1, 0 ); }

    void get_root0( T_real & re, T_real & im ) const
    {
      if ( m_cplx )
      {
        re = m_r0;
        im = m_r1;
      }
      else
      {
        re = m_r0;
        im = 0;
      }
    }

    void getRoot0( T_real & re, T_real & im ) const { return get_root0( re, im ); }

    void get_root0( T_complex & r ) const { r = m_cplx ? T_complex( m_r0, m_r1 ) : T_complex( m_r0, 0 ); }

    void getRoot0( T_complex & r ) const { return get_root0( r ); }

    void get_root1( T_real & re, T_real & im ) const
    {
      if ( m_cplx )
      {
        re = m_r0;
        im = -m_r1;
      }
      else
      {
        re = m_r1;
        im = 0;
      }
    }

    void getRoot1( T_real & re, T_real & im ) const { return get_root1( re, im ); }

    void get_root1( T_complex & r ) const { r = m_cplx ? T_complex( m_r0, -m_r1 ) : T_complex( m_r1, 0 ); }

    void getRoot1( T_complex & r ) const { return get_root1( r ); }

    T_complex root( integer const i ) const
    {
      switch ( i )
      {
        case 0: return root0();
        case 1: return root1();
      }
      return 0;
    }

    void get_root( integer const i, T_real & re, T_real & im ) const
    {
      switch ( i )
      {
        case 0: return get_root0( re, im );
        case 1: return get_root1( re, im );
      }
    }

    T_real eval( T_real const & x ) const { return eval_poly<T_real>( m_ABC, 2, x ); }

    T_complex eval( T_complex const & x ) const { return eval_poly_complex<T_real, T_complex>( m_ABC, 2, x ); }

    void eval( T_real const & x, T_real & p, T_real & dp ) const { eval_poly_Dpoly<T_real>( m_ABC, 2, x, p, dp ); }

    void info( ostream_type & s ) const;

    [[nodiscard]] bool check( ostream_type & s ) const;
  };
 
  using Quadratic   = QuadraticT<real_type, real_complex>;
#if POLYNOMIAL_ROOTS_HAS_MULTIPRECISION
  using QuadraticHQ = QuadraticT<quad_real, quad_complex>;
#endif
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  /*\
   |    ____      _     _
   |   / ___|   _| |__ (_) ___
   |  | |  | | | | '_ \| |/ __|
   |  | |__| |_| | |_) | | (__
   |   \____\__,_|_.__/|_|\___|
   |
   |  A * x^3 + B * x^2 + C * x + D
  \*/
  template <typename T_real, typename T_complex> class CubicT
  {
    T_real  m_ABCD[4]{ 0, 0, 0, 0 };
    T_real  m_r0   = 0;
    T_real  m_r1   = 0;
    T_real  m_r2   = 0;
    integer m_nrts = 0;
    integer m_iter = 0;
    bool    m_cplx = false;  // complex root
    bool    m_dblx = false;  // double root
    bool    m_trpx = false;  // triple root
 
    void find_roots();
 
  public:
    using value_type   = T_real;
    using complex_type = T_complex;

    CubicT() = default;

    CubicT( T_real const & a, T_real const & b, T_real const & c, T_real const & d )
    {
      m_ABCD[0] = a;
      m_ABCD[1] = b;
      m_ABCD[2] = c;
      m_ABCD[3] = d;
      // find roots only on finite values
      root_assert(
        isfinite( a ) && isfinite( b ) && isfinite( c ) && isfinite( d ),
        "CubicT( a={}, b={}, c={}, d={} ) arguments must be finite!",
        a, b, c, d
      );
      find_roots();
    }

    void setup( T_real const & a, T_real const & b, T_real const & c, T_real const & d )
    {
      m_ABCD[0] = a;
      m_ABCD[1] = b;
      m_ABCD[2] = c;
      m_ABCD[3] = d;
      m_nrts    = 0;
      m_iter    = 0;
      m_cplx    = false;  // complex root
      m_dblx    = false;  // double root
      m_trpx    = false;  // triple root
      root_assert(
        isfinite( a ) && isfinite( b ) && isfinite( c ) && isfinite( d ),
        "CubicT::setup( a={}, b={}, c={}, d={} ) arguments must be finite!",
        a, b, c, d
      );
      find_roots();
    }

    integer num_roots() const { return m_nrts; }

    integer numRoots() const { return m_nrts; }

    bool complex_root() const { return m_cplx; }

    bool complexRoot() const { return m_cplx; }

    bool double_root() const { return m_dblx; }

    bool doubleRoot() const { return m_dblx; }

    bool triple_root() const { return m_trpx; }

    bool tripleRoot() const { return m_trpx; }

    integer get_real_roots( T_real r[] ) const;

    integer getRealRoots( T_real r[] ) const { return get_real_roots( r ); }

    integer get_positive_roots( T_real r[] ) const;

    integer getPositiveRoots( T_real r[] ) const { return get_positive_roots( r ); }

    integer get_negative_roots( T_real r[] ) const;

    integer getNegativeRoots( T_real r[] ) const { return get_negative_roots( r ); }

    integer get_roots_in_range( T_real const & a, T_real const & b, T_real r[] ) const;

    integer getRootsInRange( T_real const & a, T_real const & b, T_real r[] ) const
    { return get_roots_in_range( a, b, r ); }

    integer get_roots_in_open_range( T_real const & a, T_real const & b, T_real r[] ) const;

    integer getRootsInOpenRange( T_real const & a, T_real const & b, T_real r[] ) const
    { return get_roots_in_open_range( a, b, r ); }

    T_real real_root0() const { return m_r0; }

    T_real real_root1() const { return m_r1; }

    T_real real_root2() const { return m_r2; }

    T_complex root0() const { return m_cplx ? T_complex( m_r0, m_r1 ) : T_complex( m_r0, 0 ); }

    T_complex root1() const { return m_cplx ? T_complex( m_r0, -m_r1 ) : T_complex( m_r1, 0 ); }

    T_complex root2() const { return T_complex( m_r2, 0 ); }

    void get_root0( T_real & re, T_real & im ) const
    {
      if ( m_cplx )
      {
        re = m_r0;
        im = m_r1;
      }
      else
      {
        re = m_r0;
        im = 0;
      }
    }

    void getRoot0( T_real & re, T_real & im ) const { get_root0( re, im ); }

    void get_root0( T_complex & r ) const { r = m_cplx ? T_complex( m_r0, m_r1 ) : T_complex( m_r0, 0 ); }

    void getRoot0( T_complex & r ) const { get_root0( r ); }

    void get_root1( T_real & re, T_real & im ) const
    {
      if ( m_cplx )
      {
        re = m_r0;
        im = -m_r1;
      }
      else
      {
        re = m_r1;
        im = 0;
      }
    }

    void getRoot1( T_real & re, T_real & im ) const { get_root1( re, im ); }

    void get_root1( T_complex & r ) const { r = m_cplx ? T_complex( m_r0, -m_r1 ) : T_complex( m_r1, 0 ); }

    void getRoot1( T_complex & r ) const { get_root1( r ); }

    void get_root2( T_real & re, T_real & im ) const
    {
      re = m_r2;
      im = 0;
    }

    void getRoot2( T_real & re, T_real & im ) const { get_root2( re, im ); }

    void get_root2( T_complex & r ) const { r = T_complex( m_r2, 0 ); }

    void getRoot2( T_complex & r ) const { get_root2( r ); }

    T_complex root( integer const i ) const
    {
      switch ( i )
      {
        case 0: return root0();
        case 1: return root1();
        case 2: return root2();
      }
      return 0;
    }

    void get_root( integer const i, T_real & re, T_real & im ) const
    {
      switch ( i )
      {
        case 0: return get_root0( re, im );
        case 1: return get_root1( re, im );
        case 2: return get_root2( re, im );
      }
    }

    T_real eval( T_real const & x ) const { return eval_poly<T_real>( m_ABCD, 3, x ); }

    T_complex eval( T_complex const & x ) const { return eval_poly_complex<T_real, T_complex>( m_ABCD, 3, x ); }

    void eval( T_real const & x, T_real & p, T_real & dp ) const { eval_poly_Dpoly<T_real>( m_ABCD, 3, x, p, dp ); }

    void info( ostream_type & s ) const;

    [[nodiscard]] bool check( ostream_type & s ) const;
  };
 
  using Cubic   = CubicT<real_type, real_complex>;
#if POLYNOMIAL_ROOTS_HAS_MULTIPRECISION
  using CubicHQ = CubicT<quad_real, quad_complex>;
#endif
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  /*\
   |    ___                   _   _
   |   / _ \ _   _  __ _ _ __| |_(_) ___
   |  | | | | | | |/ _` | '__| __| |/ __|
   |  | |_| | |_| | (_| | |  | |_| | (__
   |   \__\_\\__,_|\__,_|_|   \__|_|\___|
   |
   |  A * x^4 + B * x^3 + C * x^2 + D * x + E
  \*/
  template <typename T_real, typename T_complex> class QuarticT
  {
    T_real  m_ABCDE[5]{ 0, 0, 0, 0, 0 };
    T_real  m_r0    = 0;
    T_real  m_r1    = 0;
    T_real  m_r2    = 0;
    T_real  m_r3    = 0;
    integer m_iter  = 0;
    integer m_nreal = 0;
    integer m_ncplx = 0;
 
    void find_roots();
 
    bool cplx0() const { return m_ncplx > 0; }
    bool cplx1() const { return m_ncplx > 0; }
    bool cplx2() const { return m_ncplx > 2; }
    bool cplx3() const { return m_ncplx > 2; }
 
  public:
    using value_type   = T_real;
    using complex_type = T_complex;

    QuarticT() = default;

    QuarticT( T_real const & a, T_real const & b, T_real const & c, T_real const & d, T_real const & e )
    {
      m_ABCDE[0] = a;
      m_ABCDE[1] = b;
      m_ABCDE[2] = c;
      m_ABCDE[3] = d;
      m_ABCDE[4] = e;
      // find roots only on finite values
      root_assert(
        isfinite( a ) && isfinite( b ) && isfinite( c ) && isfinite( d ) && isfinite( e ),
        "QuarticT( a={}, b={}, c={}, d={}, e={} ) arguments must be finite!",
        a, b, c, d, e
      );
      find_roots();
    }

    void setup( T_real const & a, T_real const & b, T_real const & c, T_real const & d, T_real const & e )
    {
      m_ABCDE[0] = a;
      m_ABCDE[1] = b;
      m_ABCDE[2] = c;
      m_ABCDE[3] = d;
      m_ABCDE[4] = e;
      m_iter     = 0;
      m_nreal    = 0;
      m_ncplx    = 0;
      // find roots only on finite values
      root_assert(
        isfinite( a ) && isfinite( b ) && isfinite( c ) && isfinite( d ) && isfinite( e ),
        "QuarticT::setup( a={}, b={}, c={}, d={}, e={} ) arguments must be finite!",
        a, b, c, d, e
      );
      find_roots();
    }

    integer num_roots() const { return m_nreal + m_ncplx; }

    integer numRoots() const { return m_nreal + m_ncplx; }

    integer num_real_roots() const { return m_nreal; }

    integer numRealRoots() const { return m_nreal; }

    integer num_complex_roots() const { return m_ncplx; }

    integer numComplexRoots() const { return m_ncplx; }

    integer get_real_roots( T_real r[] ) const;

    integer getRealRoots( T_real r[] ) const { return get_real_roots( r ); }

    integer get_positive_roots( T_real r[] ) const;

    integer getPositiveRoots( T_real r[] ) const { return get_positive_roots( r ); }

    integer get_negative_roots( T_real r[] ) const;

    integer getNegativeRoots( T_real r[] ) const { return get_negative_roots( r ); }

    integer get_roots_in_range( T_real const & a, T_real const & b, T_real r[] ) const;

    integer getRootsInRange( T_real const & a, T_real const & b, T_real r[] ) const
    { return get_roots_in_range( a, b, r ); }

    integer get_roots_in_open_range( T_real const & a, T_real const & b, T_real r[] ) const;

    integer getRootsInOpenRange( T_real const & a, T_real const & b, T_real r[] ) const
    { return get_roots_in_open_range( a, b, r ); }

    T_real real_root0() const { return m_r0; }

    T_real real_root1() const { return m_r1; }

    T_real real_root2() const { return m_r2; }

    T_real real_root3() const { return m_r3; }

    T_complex root0() const { return cplx0() ? T_complex( m_r0, m_r1 ) : T_complex( m_r0, 0 ); }

    T_complex root1() const { return cplx1() ? T_complex( m_r0, -m_r1 ) : T_complex( m_r1, 0 ); }

    T_complex root2() const { return cplx2() ? T_complex( m_r2, m_r3 ) : T_complex( m_r2, 0 ); }

    T_complex root3() const { return cplx3() ? T_complex( m_r2, -m_r3 ) : T_complex( m_r3, 0 ); }

    void get_root0( T_real & re, T_real & im ) const
    {
      if ( cplx0() )
      {
        re = m_r0;
        im = m_r1;
      }
      else
      {
        re = m_r0;
        im = 0;
      }
    }

    void getRoot0( T_real & re, T_real & im ) const { get_root0( re, im ); }

    void get_root0( T_complex & r ) const
    {
      if ( cplx0() )
        r = T_complex( m_r0, m_r1 );
      else
        r = T_complex( m_r0, 0 );
    }

    void getRoot0( T_complex & r ) const { get_root0( r ); }

    void get_root1( T_real & re, T_real & im ) const
    {
      if ( cplx1() )
      {
        re = m_r0;
        im = -m_r1;
      }
      else
      {
        re = m_r1;
        im = 0;
      }
    }

    void getRoot1( T_real & re, T_real & im ) const { get_root1( re, im ); }

    void get_root1( T_complex & r ) const
    {
      if ( cplx1() )
        r = T_complex( m_r0, -m_r1 );
      else
        r = T_complex( m_r1, 0 );
    }

    void getRoot1( T_complex & r ) const { get_root1( r ); }

    void get_root2( T_real & re, T_real & im ) const
    {
      if ( cplx2() )
      {
        re = m_r2;
        im = m_r3;
      }
      else
      {
        re = m_r2;
        im = 0;
      }
    }

    void getRoot2( T_real & re, T_real & im ) const { get_root2( re, im ); }

    void get_root2( T_complex & r ) const
    {
      if ( cplx2() )
        r = T_complex( m_r2, m_r3 );
      else
        r = T_complex( m_r2, 0 );
    }

    void getRoot2( T_complex & r ) const { get_root2( r ); }

    void get_root3( T_real & re, T_real & im ) const
    {
      if ( cplx3() )
      {
        re = m_r2;
        im = -m_r3;
      }
      else
      {
        re = m_r3;
        im = 0;
      }
    }

    void getRoot3( T_real & re, T_real & im ) const { get_root3( re, im ); }

    void get_root3( T_complex & r ) const
    {
      if ( cplx3() )
        r = T_complex( m_r2, -m_r3 );
      else
        r = T_complex( m_r3, 0 );
    }

    void getRoot3( T_complex & r ) const { get_root3( r ); }

    T_complex root( integer const i ) const
    {
      switch ( i )
      {
        case 0: return root0();
        case 1: return root1();
        case 2: return root2();
        case 3: return root3();
      }
      return 0;
    }

    void get_root( integer const i, T_real & re, T_real & im ) const
    {
      switch ( i )
      {
        case 0: return get_root0( re, im );
        case 1: return get_root1( re, im );
        case 2: return get_root2( re, im );
        case 3: return get_root3( re, im );
      }
    }

    T_real eval( T_real const & x ) const { return eval_poly<T_real>( m_ABCDE, 4, x ); }

    T_complex eval( T_complex const & x ) const { return eval_poly_complex<T_real, T_complex>( m_ABCDE, 4, x ); }

    void eval( T_real const & x, T_real & p, T_real & dp ) const { eval_poly_Dpoly<T_real>( m_ABCDE, 4, x, p, dp ); }

    void info( ostream_type & s ) const;

    [[nodiscard]] bool check( ostream_type & s ) const;
  };
 
  using Quartic   = QuarticT<real_type, real_complex>;
#if POLYNOMIAL_ROOTS_HAS_MULTIPRECISION
  using QuarticHQ = QuarticT<quad_real, quad_complex>;
#endif
 
  /*\
   |   _   _ _   _ _
   |  | | | | |_(_) |___
   |  | | | | __| | / __|
   |  | |_| | |_| | \__ \
   |   \___/ \__|_|_|___/
  \*/

  template <typename T_real>
  inline T_real evalMonicCubic( T_real const & x, T_real const & a, T_real const & b, T_real const & c )
  {
    T_real p;
    p = x + a;
    p = p * x + b;
    p = p * x + c;
    return p;
  }

  template <typename T_real> inline void evalMonicCubic(
    T_real const & x,
    T_real const & a,
    T_real const & b,
    T_real const & c,
    T_real &       p,
    T_real &       dp )
  {
    p  = x + a;
    dp = x + p;
    p  = p * x + b;
    dp = dp * x + p;
    p  = p * x + c;
  }

  template <typename T_real> inline void evalMonicCubic(
    T_real const & x,
    T_real const & a,
    T_real const & b,
    T_real const & c,
    T_real &       p,
    T_real &       dp,
    T_real &       ddp )
  {
    p   = x + a;
    dp  = x + p;      // 2*x + a
    p   = p * x + b;  // x^2 + a * x + b
    ddp = 2 * ( x + dp );
    dp  = dp * x + p;
    p   = p * x + c;
  }

  template <typename T_real> inline T_real evalMonicQuartic(
    T_real const & x,
    T_real const & a,
    T_real const & b,
    T_real const & c,
    T_real const & d )
  {
    T_real p;
    p = x + a;      // x + a
    p = p * x + b;  // x^2+ a*x + b
    p = p * x + c;  // x^3+ a*x^2 + b*x + c
    p = p * x + d;  // x^4+ a*x^3 + b*x^2 + c*x + d
    return p;
  }

  template <typename T_real> inline void evalMonicQuartic(
    T_real const & x,
    T_real const & a,
    T_real const & b,
    T_real const & c,
    T_real const & d,
    T_real &       p,
    T_real &       dp )
  {
    p  = x + a;       // x + a
    dp = x + p;       // 2*x + a
    p  = p * x + b;   // x^2+ a*x + b
    dp = dp * x + p;  // 3*x^2 + 2*a*x + b
    p  = p * x + c;   // x^3+ a*x^2 + b*x + c
    dp = dp * x + p;  // 4*x^3 + 3*a*x^2 + 2*b*x + c
    p  = p * x + d;   // x^4+ a*x^3 + b*x^2 + c*x + d
  }

  template <typename T_real> inline void evalMonicQuartic(
    T_real const & x,
    T_real const & a,
    T_real const & b,
    T_real const & c,
    T_real const & d,
    T_real &       p,
    T_real &       dp,
    T_real &       ddp )
  {
    // p_{n+1}(x)   = x * p_{n}(x) + b_{n}
    // p'_{n+1}(x)  = x * p'_{n}(x) + p_{n}(x)
    // p''_{n+1}(x) = x * p''_{n}(x) + 2*p'_{n}(x)
    // ddp = 0;
    // dp  = 1;
    p = x + a;  // x + a
 
    ddp = 2;
    dp  = x + p;
    p   = p * x + b;
 
    ddp = ddp * x + 2 * dp;
    dp  = dp * x + p;
    p   = p * x + c;
 
    ddp = ddp * x + 2 * dp;
    dp  = dp * x + p;
    p   = p * x + d;
  }
 
}  // namespace PolynomialRoots
 
#endif