# Program Listing for File PolynomialRoots.hh

↰ Return to documentation for file (`PolynomialRoots.hh`)

    /*--------------------------------------------------------------------------*\
     |                                                                          |
     |  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 {

      typedef double valueType;
      typedef int    indexType;
      typedef std::complex<valueType> complexType;

      #ifndef DOXYGEN_SHOULD_SKIP_THIS

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

      static
      inline
      bool
      isInfinite( valueType x )
      { return FP_INFINITE == std::fpclassify(x); }

      static
      inline
      bool
      isNaN( valueType x )
      { return FP_NAN == std::fpclassify(x); }

      static
      inline
      bool
      isRegular( valueType x )
      { return !( FP_INFINITE == std::fpclassify(x) ||
                  FP_NAN      == std::fpclassify(x) ); }

      valueType
      evalPoly(
        valueType const op[],
        indexType       Degree,
        valueType       x
      );

      void
      evalPolyDPoly(
        valueType const op[],
        indexType       Degree,
        valueType       x,
        valueType     & p,
        valueType     & dp
      );

      bool
      NewtonStep(
        valueType const op[],
        indexType       Degree,
        valueType     & x
      );

      std::complex<valueType>
      evalPolyC(
        valueType const                 op[],
        indexType                       Degree,
        std::complex<valueType> const & x
      );

      #endif

      int
      roots(
        valueType const * op,
        indexType         Degree,
        valueType       * zeror,
        valueType       * zeroi
      );

      // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      /*\
       |    ___                  _           _   _
       |   / _ \ _   _  __ _  __| |_ __ __ _| |_(_) ___
       |  | | | | | | |/ _` |/ _` | '__/ _` | __| |/ __|
       |  | |_| | |_| | (_| | (_| | | | (_| | |_| | (__
       |   \__\_\\__,_|\__,_|\__,_|_|  \__,_|\__|_|\___|
       |
       |  A * x^2 + B * x + C
      \*/
      class Quadratic {
        valueType ABC[3];
        valueType r0, r1;
        indexType nrts;
        bool      cplx;
        bool      dblx;

        void findRoots();

      public:

        Quadratic() : nrts(0), cplx(false), dblx(false) {}

        Quadratic( valueType a, valueType b, valueType c )
        : nrts(0), cplx(false), dblx(false) {
          valueType & A = ABC[0];
          valueType & B = ABC[1];
          valueType & C = ABC[2];
          A = a; B = b; C = c;
          findRoots();
        }

        void
        setup( valueType a, valueType b, valueType c ) {
          valueType & A = ABC[0];
          valueType & B = ABC[1];
          valueType & C = ABC[2];
          A = a; B = b; C = c;
          findRoots();
        }

        indexType numRoots() const { return nrts; }

        bool complexRoots() const { return cplx; }

        bool doubleRoot() const { return dblx; }

        indexType getRealRoots( valueType r[] ) const;

        indexType getPositiveRoots( valueType r[] ) const;

        indexType getNegativeRoots( valueType r[] ) const;

        indexType getRootsInRange( valueType a, valueType b, valueType r[] ) const;

        indexType getRootsInOpenRange( valueType a, valueType b, valueType r[] ) const;

        valueType real_root0() const { return r0; }

        valueType real_root1() const { return r1; }

        complexType
        root0() const
        { return cplx ? complexType(r0,r1) : complexType(r0,0); }

        complexType
        root1() const
        { return cplx ? complexType(r0,-r1) : complexType(r1,0); }

        void
        getRoot0( valueType & re, valueType & im ) const {
          if ( cplx ) { re = r0; im = r1; }
          else        { re = r0; im = 0;  }
        }

        void
        getRoot0( complexType & r ) const
        { r = cplx ? complexType(r0,r1) : complexType(r0,0); }

        void
        getRoot1( valueType & re, valueType & im ) const {
          if ( cplx ) { re = r0; im = -r1; }
          else        { re = r1; im = 0;   }
        }

        void
        getRoot1( complexType & r ) const
        { r = cplx ? complexType(r0,-r1) : complexType(r1,0); }

        valueType
        eval( valueType x ) const
        { return evalPoly( ABC, 2, x ); }

        complexType
        eval( complexType const & x ) const
        { return evalPolyC( ABC, 2, x ); }

        void eval( valueType x, valueType & p, valueType & dp ) const;

        void
        info( std::ostream & s ) const;

        bool
        check( std::ostream & s ) const;

      };

      // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      /*\
       |    ____      _     _
       |   / ___|   _| |__ (_) ___
       |  | |  | | | | '_ \| |/ __|
       |  | |__| |_| | |_) | | (__
       |   \____\__,_|_.__/|_|\___|
       |
       |  A * x^3 + B * x^2 + C * x + D
      \*/
      class Cubic {
        valueType ABCD[4];
        valueType r0, r1, r2;
        indexType nrts, iter;
        bool      cplx; // complex root
        bool      dblx; // double root
        bool      trpx; // triple root

        void findRoots();

      public:

        Cubic() : nrts(0), iter(0), cplx(false), trpx(false) {}
        Cubic( valueType _a, valueType _b, valueType _c, valueType _d )
        : nrts(0), iter(0), cplx(false), trpx(false) {
          valueType & A = ABCD[0];
          valueType & B = ABCD[1];
          valueType & C = ABCD[2];
          valueType & D = ABCD[3];
          A = _a; B = _b; C = _c; D = _d;
          findRoots();
        }

        void
        setup( valueType a, valueType b, valueType c, valueType d ) {
          valueType & A = ABCD[0];
          valueType & B = ABCD[1];
          valueType & C = ABCD[2];
          valueType & D = ABCD[3];
          A = a; B = b; C = c; D = d;
          findRoots();
        }

        indexType numRoots()     const { return nrts; }

        bool complexRoots() const { return cplx; }

        bool doubleRoot() const { return dblx; }

        bool tripleRoot() const { return trpx; }

        indexType getRealRoots( valueType r[] ) const;

        indexType getPositiveRoots( valueType r[] ) const;

        indexType getNegativeRoots( valueType r[] ) const;

        indexType getRootsInRange( valueType a, valueType b, valueType r[] ) const;

        indexType getRootsInOpenRange( valueType a, valueType b, valueType r[] ) const;

        valueType real_root0() const { return r0; }

        valueType real_root1() const { return r1; }

        valueType real_root2() const { return r2; }

        complexType
        root0() const
        { return cplx ? complexType(r0,r1) : complexType(r0,0); }

        complexType
        root1() const
        { return cplx ? complexType(r0,-r1) : complexType(r1,0); }

        complexType
        root2() const
        { return complexType(r2,0); }

        void
        getRoot0( valueType & re, valueType & im ) const {
          if ( cplx ) { re = r0; im = r1; }
          else        { re = r0; im = 0;  }
        }

        void
        getRoot0( complexType & r ) const
        { r = cplx ? complexType(r0,r1) : complexType(r0,0); }

        void
        getRoot1( valueType & re, valueType & im ) const {
          if ( cplx ) { re = r0; im = -r1; }
          else        { re = r1; im = 0;   }
        }

        void
        getRoot1( complexType & r ) const
        { r = cplx ? complexType(r0,-r1) : complexType(r1,0); }

        void
        getRoot2( valueType & re, valueType & im ) const
        { re = r2; im = 0; }

        void
        getRoot2( complexType & r ) const
        { r = complexType(r2,0); }

        valueType
        eval( valueType x ) const
        { return evalPoly( ABCD, 3, x ); }

        complexType
        eval( complexType const & x ) const
        { return evalPolyC( ABCD, 3, x ); }

        void eval( valueType x, valueType & p, valueType & dp ) const;

        void
        info( std::ostream & s ) const;

        bool
        check( std::ostream & s ) const;

      };

      // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      /*\
       |    ___                   _   _
       |   / _ \ _   _  __ _ _ __| |_(_) ___
       |  | | | | | | |/ _` | '__| __| |/ __|
       |  | |_| | |_| | (_| | |  | |_| | (__
       |   \__\_\\__,_|\__,_|_|   \__|_|\___|
       |
       |  A * x^3 + B * x^2 + C * x + D
      \*/
      class Quartic {
        valueType ABCDE[5];
        valueType r0, r1, r2, r3;
        indexType iter, nreal, ncplx;

        void findRoots();

        bool cplx0() const { return ncplx > 0; }
        bool cplx1() const { return ncplx > 0; }
        bool cplx2() const { return ncplx > 2; }
        bool cplx3() const { return ncplx > 2; }

      public:

        Quartic() : iter(0), nreal(0), ncplx(0) {}
        Quartic(
          valueType _a,
          valueType _b,
          valueType _c,
          valueType _d,
          valueType _e
        )
        : iter(0), nreal(0), ncplx(0) {
          valueType & A = ABCDE[0];
          valueType & B = ABCDE[1];
          valueType & C = ABCDE[2];
          valueType & D = ABCDE[3];
          valueType & E = ABCDE[4];
          A = _a; B = _b; C = _c; D = _d; E = _e;
          findRoots();
        }

        void
        setup(
          valueType a,
          valueType b,
          valueType c,
          valueType d,
          valueType e
        ) {
          valueType & A = ABCDE[0];
          valueType & B = ABCDE[1];
          valueType & C = ABCDE[2];
          valueType & D = ABCDE[3];
          valueType & E = ABCDE[4];
          A = a; B = b; C = c; D = d; E = e;
          findRoots();
        }

        indexType numRoots() const { return nreal+ncplx; }

        indexType numRealRoots() const { return nreal; }

        indexType numComplexRoots() const { return ncplx; }

        indexType getRealRoots( valueType r[] ) const;

        indexType getPositiveRoots( valueType r[] ) const;

        indexType getNegativeRoots( valueType r[] ) const;

        indexType getRootsInRange( valueType a, valueType b, valueType r[] ) const;

        indexType getRootsInOpenRange( valueType a, valueType b, valueType r[] ) const;

        valueType real_root0() const { return r0; }

        valueType real_root1() const { return r1; }

        valueType real_root2() const { return r2; }

        valueType real_root3() const { return r3; }

        complexType
        root0() const
        { return cplx0() ? complexType(r0,r1) : complexType(r0,0); }

        complexType
        root1() const
        { return cplx1() ? complexType(r0,-r1) : complexType(r1,0); }

        complexType
        root2() const
        { return cplx2() ? complexType(r2,r3) : complexType(r2,0); }

        complexType
        root3() const
        { return cplx3() ? complexType(r2,-r3) : complexType(r3,0); }

        void
        getRoot0( valueType & re, valueType & im ) const {
          if ( cplx0() ) { re = r0; im = r1; }
          else           { re = r0; im = 0;  }
        }

        void
        getRoot0( complexType & r ) const {
          if ( cplx0() ) r = complexType(r0,r1);
          else           r = complexType(r0,0);
        }

        void
        getRoot1( valueType & re, valueType & im ) const {
          if ( cplx1() ) { re = r0; im = -r1; }
          else           { re = r1; im = 0;   }
        }

        void
        getRoot1( complexType & r ) const {
          if ( cplx1() ) r = complexType(r0,-r1);
          else           r = complexType(r1,0);
        }

        void
        getRoot2( valueType & re, valueType & im ) const {
          if ( cplx2() ) { re = r2; im = r3; }
          else           { re = r2; im = 0;  }
        }

        void
        getRoot2( complexType & r ) const {
          if ( cplx2() ) r = complexType(r2,r3);
          else           r = complexType(r2,0);
        }

        void
        getRoot3( valueType & re, valueType & im ) const {
          if ( cplx3() ) { re = r2; im = -r3; }
          else           { re = r3; im = 0;   }
        }

        void
        getRoot3( complexType & r ) const {
          if ( cplx3() ) r = complexType(r2,-r3);
          else           r = complexType(r3,0);
        }

        valueType
        eval( valueType x ) const
        { return evalPoly( ABCDE, 4, x ); }

        complexType
        eval( complexType const & x ) const
        { return evalPolyC( ABCDE, 4, x ); }

        void
        info( std::ostream & s ) const;

        bool
        check( std::ostream & s ) const;

      };

      /*\
       |   _   _ _   _ _
       |  | | | | |_(_) |___
       |  | | | | __| | / __|
       |  | |_| | |_| | \__ \
       |   \___/ \__|_|_|___/
      \*/

      #ifndef DOXYGEN_SHOULD_SKIP_THIS

      // x^3 + a*x^2 + b*x + c
      static
      inline
      valueType
      evalMonicCubic(
        valueType x,
        valueType a,
        valueType b,
        valueType c
      ) {
        valueType p;
        p = x + a;
        p = p * x + b;
        p = p * x + c;
        return p;
      }

      static
      inline
      void
      evalMonicCubic(
        valueType   x,
        valueType   a,
        valueType   b,
        valueType   c,
        valueType & p,
        valueType & 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(
        valueType   x,
        valueType   a,
        valueType   b,
        valueType   c,
        valueType & p,
        valueType & dp,
        valueType & 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
      valueType
      evalMonicQuartic(
        valueType x,
        valueType a,
        valueType b,
        valueType c,
        valueType d
      ) {
        valueType 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(
        valueType   x,
        valueType   a,
        valueType   b,
        valueType   c,
        valueType   d,
        valueType & p,
        valueType & 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(
        valueType   x,
        valueType   a,
        valueType   b,
        valueType   c,
        valueType   d,
        valueType & p,
        valueType & dp,
        valueType & 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;
      }

      #endif

    }

    #endif