Program Listing for File PolynomialRoots-1-Quadratic.cc

↰ Return to documentation for file (PolynomialRoots-1-Quadratic.cc)

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2017                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Universita` degli Studi di Trento                                   |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/

#ifdef __GNUC__
#pragma GCC diagnostic ignored "-Wunused-function"
#endif
#ifdef __clang__
#pragma clang diagnostic ignored "-Wglobal-constructors"
#pragma clang diagnostic ignored "-Wvla-extension"
#pragma clang diagnostic ignored "-Wvla"
#pragma clang diagnostic ignored "-Wunused-function"
#pragma clang diagnostic ignored "-Wdocumentation-unknown-command"
#endif

#include "PolynomialRoots.hh"

#include <cmath>
#include <iostream>
#include <algorithm>
#include <limits>

#ifndef DOXYGEN_SHOULD_SKIP_THIS

namespace PolynomialRoots {

  using  std::abs;
  static valueType const machepsi = std::numeric_limits<valueType>::epsilon();

  indexType
  Quadratic::getRealRoots( valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx ) {
      r[nr++] = r0;
      if ( nrts > 1 ) r[nr++] = r1;
    }
    return nr;
  }

  indexType
  Quadratic::getPositiveRoots( valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx ) {
      if ( nrts > 0 && r0 > 0 ) r[nr++] = r0;
      if ( nrts > 1 && r1 > 0 ) r[nr++] = r1;
    }
    return nr;
  }

  indexType
  Quadratic::getNegativeRoots( valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx ) {
      if ( nrts > 0 && r0 < 0 ) r[nr++] = r0;
      if ( nrts > 1 && r1 < 0 ) r[nr++] = r1;
    }
    return nr;
  }

  indexType
  Quadratic::getRootsInRange( valueType a, valueType b, valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx ) {
      if ( nrts > 0 && r0 >= a && r0 <= b ) r[nr++] = r0;
      if ( nrts > 1 && r1 >= a && r1 <= b ) r[nr++] = r1;
    }
    return nr;
  }

  indexType
  Quadratic::getRootsInOpenRange( valueType a, valueType b, valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx ) {
      if ( nrts > 0 && r0 > a && r0 < b ) r[nr++] = r0;
      if ( nrts > 1 && r1 > a && r1 < b ) r[nr++] = r1;
    }
    return nr;
  }

  void
  Quadratic::eval( valueType x, valueType & p, valueType & dp ) const {
    valueType const & A = ABC[0];
    valueType const & B = ABC[1];
    valueType const & C = ABC[2];
    if ( std::abs(x) > 1 ) {
      valueType z  = 1/x;
      valueType x2 = x*x;
      p = ((C*z+B)*z+A)*x2;
    } else {
      p = (A*x+B)*x+C;
    }
    dp = 2*A*x+B;
  }

  /*\
   *  Calculate the zeros of the quadratic a*z^2 + b*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.
  \*/

  void
  Quadratic::findRoots() {
    valueType const & A = ABC[0];
    valueType const & B = ABC[1];
    valueType const & C = ABC[2];

    r0 = r1 = 0;
    nrts = 0;
    cplx = dblx = false;

    if ( isZero(A) ) { // less than two roots b*z + c = 0
      if ( !isZero(B) ) { nrts = 1; r0 = -C/B; }
    } else if ( isZero(C) ) { // a*z^2 + b*z  = 0
      nrts = 2;
      dblx = isZero(B);
      if ( !dblx ) {
        r0 = -B/A;
        if ( r0 < 0 ) std::swap(r0,r1);
      }
    } else { // Compute discriminant avoiding overflow.
      valueType hb    = B/2; // b now b/2
      valueType abs_b = std::abs(hb);
      valueType abs_c = std::abs(C);
      valueType e, d;
      if ( abs_b < abs_c ) {
        e = C < 0 ? -A : A;
        e = (hb*hb)-e*abs_c;
        d = std::sqrt(std::abs(e));
      } else {
        e = 1 - (A/hb)*(C/hb);
        d = std::sqrt(std::abs(e))*abs_b;
      }
      nrts = 2;
      cplx = e < 0;
      if ( cplx ) {
        // complex conjugate zeros
        r0 = -hb/A;         // real part
        r1 = std::abs(d/A); // immaginary part
      } else {
        // real zeros
        dblx = isZero(d);
        if ( dblx ) {
          r0 = r1 = -hb/A;
        } else {
          if ( hb >= 0 ) d = -d;
          r0 = (d-hb)/A;
          //r1 = (-d-hb)/a;
          if ( !isZero(r0) ) r1 = (C/r0)/A;
          if ( r0 > r1 ) std::swap(r0,r1); // order roots
        }
      }
    }
  }

  void
  Quadratic::info( std::ostream & s ) const {
    valueType const & A = ABC[0];
    valueType const & B = ABC[1];
    valueType const & C = ABC[2];
    s << "\npoly A=" << A << " B=" << B << " C=" << C
      << "\nn. roots = " << nrts
      << "\ncomplex  = " << (cplx?"YES":"NO")
      << "\ndouble   = " << (dblx?"YES":"NO");
    if ( cplx ) {
      s << "\nx0 = (" << r0 << "," <<  r1 << ')'
        << "\nx1 = (" << r0 << "," << -r1 << ')';
    } else if ( dblx ) {
      s << "\nx0 = x1 = " << r0;
    } else if ( nrts == 1 ) {
      s << "\nx0 = " << r0;
    } else if ( nrts == 2 ) {
      s << "\nx0 = " << r0
        << "\nx1 = " << r1;
    }
    s << '\n';
  }

  bool
  Quadratic::check( std::ostream & s ) const {
    valueType const & A = ABC[0];
    valueType const & B = ABC[1];
    valueType const & C = ABC[2];
    bool ok = true;
    valueType epsi = 10 * ( std::abs(A) +
                            std::abs(B) +
                            std::abs(C) ) * machepsi;
    if ( cplx ) {
      valueType z0 = std::abs(eval( root0() ));
      valueType z1 = std::abs(eval( root1() ));
      s << "|p(r0)| = " << z0
        << "\n|p(r1)| = " << z1
        << '\n';
      ok = z0 < epsi && z1 < epsi;
    } else if ( nrts == 1 ) {
      valueType z0 = eval( real_root0() );
      s << "p(r0) = " << z0  << '\n';
      ok = std::abs(z0) < epsi;
    } else if ( nrts == 2 ) {
      valueType z0 = eval( real_root0() );
      valueType z1 = eval( real_root1() );
      s << "p(r0) = " << z0
        << "\np(r1) = " << z1
        << '\n';
      ok = std::abs(z0) < epsi && std::abs(z1) < epsi;
    }
    return ok;
  }

}

#endif

// EOF: PolynomialRoots-1-Quadratic.cc