PolynomialRoots.hh

Go to the documentation of this file.

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2014                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Universita` degli Studi di Trento                                   |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/
 
#ifndef POLYNOMIAL_ROOTS_HH
#define POLYNOMIAL_ROOTS_HH
 
#include <cmath>
#include <cfloat>
#include <complex>
#include <iostream>

namespace PolynomialRoots {
 
  using real_type    = double;
  using integer      = int;
  using complex_type = std::complex<real_type>;
  using ostream_type = std::basic_ostream<char>;
  using istream_type = std::basic_istream<char>;
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS

  static
  inline
  bool
  isZero( real_type x )
  { return FP_ZERO == std::fpclassify(x); }

  real_type
  evalPoly(
    real_type const op[],
    integer         Degree,
    real_type       x
  );
 
  void
  evalPolyDPoly(
    real_type const op[],
    integer         Degree,
    real_type       x,
    real_type     & p,
    real_type     & dp
  );
 
  bool
  NewtonStep(
    real_type const op[],
    integer         Degree,
    real_type     & x
  );

  complex_type
  evalPolyC(
    real_type const op[],
    integer         Degree,
    complex_type    x
  );
 
  #endif

  int
  roots(
    real_type const * op,
    integer           Degree,
    real_type       * zeror,
    real_type       * zeroi
  );
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  /*\
   |    ___                  _           _   _
   |   / _ \ _   _  __ _  __| |_ __ __ _| |_(_) ___
   |  | | | | | | |/ _` |/ _` | '__/ _` | __| |/ __|
   |  | |_| | |_| | (_| | (_| | | | (_| | |_| | (__
   |   \__\_\\__,_|\__,_|\__,_|_|  \__,_|\__|_|\___|
   |
   |  A * x^2 + B * x + C
  \*/
  class Quadratic {
    real_type m_ABC[3]{0,0,0};
    real_type m_r0{0}, m_r1{0};
    integer   m_nrts{0};
    bool      m_cplx{false};
    bool      m_dblx{false};
 
    void find_roots();
 
  public:
 
    Quadratic() {}

    Quadratic( real_type a, real_type b, real_type c ) {
      using std::isfinite;
      real_type & A{ m_ABC[0] };
      real_type & B{ m_ABC[1] };
      real_type & C{ m_ABC[2] };
      A = a; B = b; C = c;
      // find roots only on finite values
      if ( isfinite(a) && isfinite(b) && isfinite(c) ) {
        find_roots();
      }
    }

    void
    setup( real_type a, real_type b, real_type c ) {
      real_type & A = m_ABC[0];
      real_type & B = m_ABC[1];
      real_type & C = m_ABC[2];
      A = a; B = b; C = 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( real_type r[] ) const;

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

    integer get_positive_roots( real_type r[] ) const;

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

    integer get_negative_roots( real_type r[] ) const;

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

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

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

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

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

    real_type real_root0() const { return m_r0; }

    real_type real_root1() const { return m_r1; }

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

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

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

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

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

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

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

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

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

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

    real_type
    eval( real_type x ) const
    { return evalPoly( m_ABC, 2, x ); }

    complex_type
    eval( complex_type x ) const
    { return evalPolyC( m_ABC, 2, x ); }

    void
    eval( real_type x, real_type & p, real_type & dp ) const {
      evalPolyDPoly( m_ABC, 2, x, p, dp );
    }

    void
    info( ostream_type & s ) const;

    bool
    check( ostream_type & s ) const;
  };
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  /*\
   |    ____      _     _
   |   / ___|   _| |__ (_) ___
   |  | |  | | | | '_ \| |/ __|
   |  | |__| |_| | |_) | | (__
   |   \____\__,_|_.__/|_|\___|
   |
   |  A * x^3 + B * x^2 + C * x + D
  \*/
  class Cubic {
    real_type m_ABCD[4]{0,0,0,0};
    real_type m_r0{0}, m_r1{0}, 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:

    Cubic() {}

    Cubic( real_type a, real_type b, real_type c, real_type d ) {
      using std::isfinite;
      real_type & A = m_ABCD[0];
      real_type & B = m_ABCD[1];
      real_type & C = m_ABCD[2];
      real_type & D = m_ABCD[3];
      A = a; B = b; C = c; D = d;
      // find roots only on finite values
      if ( isfinite(a) && isfinite(b) && isfinite(c) && isfinite(d) ) {
        find_roots();
      }
    }

    void
    setup( real_type a, real_type b, real_type c, real_type d ) {
      using std::isfinite;
      real_type & A{m_ABCD[0]};
      real_type & B{m_ABCD[1]};
      real_type & C{m_ABCD[2]};
      real_type & D{m_ABCD[3]};
      m_nrts = 0;
      m_iter = 0;
      m_cplx = false; // complex root
      m_dblx = false; // double root
      m_trpx = false; // triple root
      A = a; B = b; C = c; D = d;
      if ( isfinite(a) && isfinite(b) && isfinite(c) && isfinite(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( real_type r[] ) const;

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

    integer get_positive_roots( real_type r[] ) const;

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

    integer get_negative_roots( real_type r[] ) const;

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

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

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

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

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

    real_type real_root0() const { return m_r0; }

    real_type real_root1() const { return m_r1; }

    real_type real_root2() const { return m_r2; }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    real_type
    eval( real_type x ) const
    { return evalPoly( m_ABCD, 3, x ); }

    complex_type
    eval( complex_type x ) const
    { return evalPolyC( m_ABCD, 3, x ); }

    void
    eval( real_type x, real_type & p, real_type & dp ) const {
      evalPolyDPoly( m_ABCD, 3, x, p, dp );
    }

    void
    info( ostream_type & s ) const;

    bool
    check( ostream_type & s ) const;
  };
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  /*\
   |    ___                   _   _
   |   / _ \ _   _  __ _ _ __| |_(_) ___
   |  | | | | | | |/ _` | '__| __| |/ __|
   |  | |_| | |_| | (_| | |  | |_| | (__
   |   \__\_\\__,_|\__,_|_|   \__|_|\___|
   |
   |  A * x^4 + B * x^3 + C * x^2 + D * x + E
  \*/
  class Quartic {
    real_type m_ABCDE[5]{0,0,0,0,0};
    real_type m_r0{0}, m_r1{0}, m_r2{0},m_r3{0};
    integer   m_iter{0}, m_nreal{0}, 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:
 
    Quartic() {}
    Quartic(
      real_type a,
      real_type b,
      real_type c,
      real_type d,
      real_type e
    ) {
      using std::isfinite;
      real_type & A{m_ABCDE[0]};
      real_type & B{m_ABCDE[1]};
      real_type & C{m_ABCDE[2]};
      real_type & D{m_ABCDE[3]};
      real_type & E{m_ABCDE[4]};
      A = a; B = b; C = c; D = d; E = e;
      // find roots only on finite values
      if ( isfinite(a) && isfinite(b) && isfinite(c) && isfinite(d) && isfinite(e) ) {
        find_roots();
      }
    }

    void
    setup(
      real_type a,
      real_type b,
      real_type c,
      real_type d,
      real_type e
    ) {
      using std::isfinite;
      real_type & A{m_ABCDE[0]};
      real_type & B{m_ABCDE[1]};
      real_type & C{m_ABCDE[2]};
      real_type & D{m_ABCDE[3]};
      real_type & E{m_ABCDE[4]};
      A = a; B = b; C = c; D = d; E = e;
      m_iter  = 0;
      m_nreal = 0;
      m_ncplx = 0;
      // find roots only on finite values
      if ( isfinite(a) && isfinite(b) && isfinite(c) && isfinite(d) && isfinite(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( real_type r[] ) const;

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

    integer get_positive_roots( real_type r[] ) const;

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

    integer get_negative_roots( real_type r[] ) const;

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

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

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

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

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

    real_type real_root0() const { return m_r0; }

    real_type real_root1() const { return m_r1; }

    real_type real_root2() const { return m_r2; }

    real_type real_root3() const { return m_r3; }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    real_type
    eval( real_type x ) const
    { return evalPoly( m_ABCDE, 4, x ); }

    complex_type
    eval( complex_type x ) const
    { return evalPolyC( m_ABCDE, 4, x ); }

    void
    eval( real_type x, real_type & p, real_type & dp ) const {
      evalPolyDPoly( m_ABCDE, 4, x, p, dp );
    }

    void
    info( ostream_type & s ) const;

    bool
    check( ostream_type & s ) const;
 
  };
 
  /*\
   |   _   _ _   _ _
   |  | | | | |_(_) |___
   |  | | | | __| | / __|
   |  | |_| | |_| | \__ \
   |   \___/ \__|_|_|___/
  \*/
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS
 
  // x^3 + a*x^2 + b*x + c
  static
  inline
  real_type
  evalMonicCubic(
    real_type x,
    real_type a,
    real_type b,
    real_type c
  ) {
    real_type p;
    p = x + a;
    p = p * x + b;
    p = p * x + c;
    return p;
  }
 
  static
  inline
  void
  evalMonicCubic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type & p,
    real_type & dp
  ) {
    p  = x + a;
    dp = x + p;
    p  = p  * x + b;
    dp = dp * x + p;
    p  = p  * x + c;
  }
 
  // 3*x^2 + 2*a*x + b
  // 6*x + 2*a
  static
  inline
  void
  evalMonicCubic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type & p,
    real_type & dp,
    real_type & 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;
  }
 
  // x^4 + a*x^3 + b*x^2 + c*x + d
  static
  inline
  real_type
  evalMonicQuartic(
    real_type x,
    real_type a,
    real_type b,
    real_type c,
    real_type d
  ) {
    real_type 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;
  }
 
  static
  inline
  void
  evalMonicQuartic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type   d,
    real_type & p,
    real_type & 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
  }
 
  static
  inline
  void
  evalMonicQuartic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type   d,
    real_type & p,
    real_type & dp,
    real_type & 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;
  }
 
  // x^3 + a*x^2 + b*x + c
  static
  inline
  real_type
  eval_monic_cubic(
    real_type x,
    real_type a,
    real_type b,
    real_type c
  ) {
    return evalMonicCubic( x, a, b, c );
  }
 
  static
  inline
  void
  eval_monic_cubic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type & p,
    real_type & dp
  ) {
    evalMonicCubic( x, a, b, c, p, dp );
  }
 
  // 3*x^2 + 2*a*x + b
  // 6*x + 2*a
  static
  inline
  void
  eval_monic_cubic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type & p,
    real_type & dp,
    real_type & ddp
  ) {
    evalMonicCubic( x, a, b, c, p, dp, ddp );
  }
 
  // x^4 + a*x^3 + b*x^2 + c*x + d
  static
  inline
  real_type
  eval_monic_quartic(
    real_type x,
    real_type a,
    real_type b,
    real_type c,
    real_type d
  ) {
    return evalMonicQuartic( x, a, b, c, d );
  }
 
  static
  inline
  void
  eval_monic_quartic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type   d,
    real_type & p,
    real_type & dp
  ) {
    evalMonicQuartic( x, a, b, c, d, p, dp );
  }
 
  static
  inline
  void
  eval_monic_quartic(
    real_type   x,
    real_type   a,
    real_type   b,
    real_type   c,
    real_type   d,
    real_type & p,
    real_type & dp,
    real_type & ddp
  ) {
    evalMonicQuartic( x, a, b, c, d, p, dp, ddp );
  }
  #endif
}
 
#endif