QGpCoreTools/RootSolver.h

Go to the documentation of this file.

/***************************************************************************
**
**  This file is part of QGpCoreTools.
**
**  This library is free software; you can redistribute it and/or
**  modify it under the terms of the GNU Lesser General Public
**  License as published by the Free Software Foundation; either
**  version 2.1 of the License, or (at your option) any later version.
**
**  This file is distributed in the hope that it will be useful, but WITHOUT
**  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
**  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
**  License for more details.
**
**  You should have received a copy of the GNU Lesser General Public
**  License along with this library; if not, write to the Free Software
**  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
**
**  See http://www.geopsy.org for more information.
**
**  Created : 2004-04-28
**  Authors :
**    Marc Wathelet
**    Marc Wathelet (ULg, Liège, Belgium)
**    Marc Wathelet (LGIT, Grenoble, France)
**
***************************************************************************/

#ifndef ROOTSOLVER_H
#define ROOTSOLVER_H

#include "QGpCoreToolsDLLExport.h"
#include "CoreApplication.h"

namespace QGpCoreTools {

#define ROOTSOLVER_MAX_ITERATIONS 500

template <class RealType, class FunctionClass>
class QGPCORETOOLS_EXPORT RootSolverTemplate
{
public:
  RootSolverTemplate(FunctionClass * values);

  FunctionClass * values() {return _values;}

  RealType searchStep() const {return _dx;}
  void setAbsoluteStep(RealType step) {_dx=fabs(step); _dxType=Absolute;}
  void setRelativeStep(RealType step) {_dx=fabs(step); _dxType=Relative;}

  void setPolarity(RealType x0);
  void inversePolarity() {_polarity=-_polarity;}

  RealType precision() const {return _precision;}
  void setPrecision(RealType prec) {_precision=prec;}

  inline bool searchUp(RealType from, RealType to);
  inline bool searchDown(RealType from, RealType to);

  inline bool searchUp(RealType from, RealType min, RealType max);
  inline bool searchDown(RealType from, RealType min, RealType max);

  inline bool searchUpSlope(RealType from, RealType min, RealType max);
  inline bool searchDownSlope(RealType from, RealType min, RealType max);

  bool setInterval(RealType min, RealType max);

  void halving();
  void newtonRaphson();
  void neville();

  RealType lower() const {return _x1<_x2 ? _x1 : _x2;}
  RealType upper() const {return _x1<_x2 ? _x2 : _x1;}
private:
  FunctionClass * _values;
  RealType _x1, _x2, _y1, _y2;

  enum StepType {Absolute, Relative};
  StepType _dxType;
  RealType _dx;
  RealType _polarity;
  RealType _precision;

  template <class Iterator>
    bool search(RealType from, RealType to, const Iterator& it);
  template <class Iterator>
    bool search(RealType from, RealType min, RealType max, const Iterator& it);
  template <class Iterator>
    bool searchSlope(RealType from, RealType min, RealType max, const Iterator& it);
  static inline bool isFound(RealType y1, RealType y2);

  class AbsoluteUpIterator
  {
  public:
    AbsoluteUpIterator(RealType dx) {_dx=dx;}
    RealType next(RealType x) const {return x+_dx;}
    static bool atEnd(RealType x, RealType to) {return x >= to;}
    static bool atEnd(RealType x, RealType /*min*/, RealType max) {return x >= max;}
    static RealType begin(RealType min, RealType /*max*/ ) {return min;}
    static RealType end(RealType /*min*/, RealType max) {return max;}
  private:
    RealType _dx;
  };

  class AbsoluteDownIterator
  {
  public:
    AbsoluteDownIterator(RealType dx) {_dx=dx;}
    RealType next(RealType x) const {return x-_dx;}
    static bool atEnd(RealType x, RealType to) {return x <= to;}
    static bool atEnd(RealType x, RealType min, RealType /*max*/ ) {return x <= min;}
    static RealType begin(RealType /*min*/, RealType max) {return max;}
    static RealType end(RealType min, RealType /*max*/ ) {return min;}
  private:
    RealType _dx;
  };

  class RelativeUpIterator
  {
  public:
    RelativeUpIterator(RealType dx) {_dx=1.0+dx;}
    RealType next(RealType x) const {return x*_dx;}
    static bool atEnd(RealType x, RealType to) {return x >= to;}
    static bool atEnd(RealType x, RealType /*min*/, RealType max) {return x >= max;}
    static RealType begin(RealType min, RealType /*max*/ ) {return min;}
    static RealType end(RealType /*min*/, RealType max) {return max;}
  private:
    RealType _dx;
  };

  class RelativeDownIterator
  {
  public:
    RelativeDownIterator(RealType dx) {_dx=1.0-dx;}
    RealType next(RealType x) const {return x*_dx;}
    static bool atEnd(RealType x, RealType to) {return x <= to;}
    static bool atEnd(RealType x, RealType min, RealType /*max*/ ) {return x <= min;}
    static RealType begin(RealType /*min*/, RealType max) {return max;}
    static RealType end(RealType min, RealType /*max*/ ) {return min;}
  private:
    RealType _dx;
  };
};

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>
::searchDown(RealType from, RealType min, RealType max)
{
  if(_dxType==Absolute) {
    AbsoluteDownIterator it(_dx);
    return search(from, min, max, it);
  } else {
    RelativeDownIterator it(_dx);
    return search(from, min, max, it);
  }
}

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>
::searchUp(RealType from, RealType min, RealType max)
{
  if(_dxType==Absolute) {
    AbsoluteUpIterator it(_dx);
    return search(from, min, max, it);
  } else {
    RelativeUpIterator it(_dx);
    return search(from, min, max, it);
  }
}

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>
::searchDownSlope(RealType from, RealType min, RealType max)
{
  if(_dxType==Absolute) {
    AbsoluteDownIterator it(_dx);
    return searchSlope(from, min, max, it);
  } else {
    RelativeDownIterator it(_dx);
    return searchSlope(from, min, max, it);
  }
}

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>
::searchUpSlope(RealType from, RealType min, RealType max)
{
  if(_dxType==Absolute) {
    AbsoluteUpIterator it(_dx);
    return searchSlope(from, min, max, it);
  } else {
    RelativeUpIterator it(_dx);
    return searchSlope(from, min, max, it);
  }
}

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>
::searchUp(RealType from, RealType to)
{
  if(from>to) return false;
  if(_dxType==Absolute) {
    AbsoluteUpIterator it(_dx);
    return search(from, to, it);
  } else {
    RelativeUpIterator it(_dx);
    return search(from, to, it);
  }
}

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>
::searchDown(RealType from, RealType to)
{
  if(from>to) return false;
  if(_dxType==Absolute) {
    AbsoluteDownIterator it(_dx);
    return search(from, to, it);
  } else {
    RelativeDownIterator it(_dx);
    return search(from, to, it);
  }
}

template <class RealType, class FunctionClass>
inline bool RootSolverTemplate<RealType, FunctionClass>::isFound(RealType y1, RealType y2)
{
  return (y1<0.0 && y2>0.0) || (y1>0.0 && y2<0.0);
}

// Usual precision
template <class FunctionClass>
class QGPCORETOOLS_EXPORT RootSolver: public RootSolverTemplate<double, FunctionClass>
{
public:
  RootSolver(FunctionClass * values) :
    RootSolverTemplate<double, FunctionClass>(values) {}
};

#ifdef BIGNUM
// Arbitrary precision for special cases or experimental purposes
template <class FunctionClass>
class QGPCORETOOLS_EXPORT RootSolverBigNum: public RootSolverTemplate<mp_real, FunctionClass>
{
public:
  RootSolverBigNum(FunctionClass * values) :
    RootSolverTemplate<mp_real, FunctionClass>(values) {}
};
#endif

template <class RealType, class FunctionClass>
RootSolverTemplate<RealType, FunctionClass>::
RootSolverTemplate(FunctionClass * values)
{
  _values=values;
  _dxType=Relative;
  _dx=0.1;
  _polarity=1.0;
  _precision=1e-7;
}

template <class RealType, class FunctionClass>
void RootSolverTemplate<RealType, FunctionClass>::
setPolarity(RealType x0)
{
  _polarity=_values->y(x0);
}

template <class RealType, class FunctionClass> template <class Iterator>
bool RootSolverTemplate<RealType, FunctionClass>::
search(RealType from, RealType to, const Iterator& it)
{
  //#define searchFromTo
#ifdef searchFromTo
  printf( "RootSolverTemplate::searchFromTo\n" );
#endif
  _x1=from;
  _y1=_values->y(_x1);
  _x2=it.next(_x1);
  while( ! it.atEnd(_x2, to)) {
    _y2=_values->y(_x2);
#ifdef searchFromTo
    printSearchDebug();
#endif
    if(isFound(_y1, _y2)) {
      return true;
    } else {
      _x1=_x2;
      _y1=_y2;
      _x2=it.next(_x1);
    }
  }

  // reaching the maximum slowness, tests if there's a last root
#ifdef searchFromTo
  printSearchDebug();
#endif
  if( ! it.atEnd(_x1, to)) {
    _x2=to;
    _y2=_values->y(_x2);
    if(isFound(_y1, _y2)) {
      return true;
    }
  }
  return false;
}

template <class RealType, class FunctionClass> template <class Iterator>
bool RootSolverTemplate<RealType, FunctionClass>::
search(RealType from, RealType min, RealType max, const Iterator& it)
{
  //#define searchMinMax
#ifdef searchMinMax
  printf( "RootSolverTemplate::searchMinMax\n" );
  printf( "Current polarity %lg\n", ::toDouble(_polarity));
#endif
  // Find the search direction by looking at polarity
  _x1=from;
  _y1=_values->y(_x1);
  _y2=_polarity;

  if(isFound(_y1, _y2)) {
    // Root is situated between tail and max/min (down/up)
#ifdef searchMinMax
    //printf("1/_x1=%.8lf _y1=%lg 1/_x2=%.8lf _y2=%lg\n",1/_x1,_y1,1/_x2,_y2);
    printf( "Root above _x1=%.8lf _y1=%lg _x2=%.8lf _y2=%lg\n", _x1, _y1, _x2, _y2);
#endif
    if(it.atEnd(from, min, max)) {
      /* If _x1 reaches the min (there is necessarly a solution
         greater than min because signs are different), there is a lot of
         chance that a mode jumping occured. The only one solution is to
         reduce the step size and recalculating all the curve until not
         getting this condition. */
      App::stream() << "** Warning ** : directly at the end" << endl;
      return false;
    } else {
      _x2=_x1;
      _y2=_y1;
      _x1=it.begin(min, max);
      _y1=_polarity;
      return true;
    }
  } else {
    // Root is situated below x0, between _x1 and min, thus search by step
    _x2=it.next(_x1);
    while( ! it.atEnd(_x2, min, max)) {
      _y2=_values->y(_x2);
#ifdef searchMinMax
      printSearchDebug()
#endif
      if(isFound(_y1, _y2)) {
        return true;
      } else {
        _x1=_x2;
        _y1=_y2;
        _x2=it.next(_x1);
      }
    }
  }
  // reaching the minimum slowness, tests if there's a last root
  if( ! it.atEnd(_x1, min, max)) {
    _x2=it.end(min, max);
    _y2=_values->y(_x2);
    if(isFound(_y1, _y2)) {
      return true;
    }
  }
#ifdef searchMinMax
  printf( "*** no root\n" );
#endif
  // If x1 reaches the min/max (there is maybe a solution
  // greater/less than min/max, signs are equal), there is a lot of
  // chance that a mode jumping occured for fundamental mode. The only one solution is to
  // reduce the step size and recalculating all the curve until not
  // getting this condition. For higher modes this condition can occur.
  return false;
}

template <class RealType, class FunctionClass> template <class Iterator>
bool RootSolverTemplate<RealType, FunctionClass>
::searchSlope(RealType from, RealType min, RealType max, const Iterator& it)
{
  _x1=from;
  _y1=_values->y(_x1);
  _x2=it.next(_x1);
  RealType maxDiff;
  if(_dxType==Absolute) {
    maxDiff=1e-7*_dx;
  } else {
    maxDiff=0.5e-7*fabs(_x1 + _x2) * _dx;
  }
  while( ! it.atEnd(_x2, min, max)) {
    _y2=_values->y(_x2);
    if(isFound(_y1, _y2)) {
      return true;
    }
    // A third point to get the slopes
    RealType x3=0.5 * (_x1 + _x2);
    RealType y3=_values->y(x3);
    if(y3 < _y1 && y3 < _y2) {
      // Refine the minimum
      RealType x1=_x1;
      RealType y1=_y1;
      RealType x2=_x2;
      RealType y2=_y2;
      RealType x13, y13, x32, y32;
      while(fabs(x2 - x1) > maxDiff) {
        x13=0.5 * (x1 + x3);
        y13=_values->y(x13);
        x32=0.5 * (x3 + x2);
        y32=_values->y(x32);
        if(y3 < y13 && y3 < y32) {
          x1=x13;
          y1=y13;
          x2=x32;
          y2=y32;
        } else if(y13 < y32) {
          x2=x3;
          y2=y3;
          x3=x13;
          y3=y13;
        } else {
          x1=x13;
          y1=y13;
          x3=x32;
          y3=y32;
        }
      }
      if(isFound(y3, _y2)) {
        // two hidden roots are found, return the first one
        _x2=x3;
        _y2=y3;
        return true;
      }
      // minimum without root restore and go further
    } else if(y3 > _y1 && y3 >  _y2) {
      // Refine the maximum
      RealType x1=_x1;
      RealType y1=_y1;
      RealType x2=_x2;
      RealType y2=_y2;
      RealType x13, y13, x32, y32;
      while(fabs(x1 - x2) > maxDiff) {
        x13=0.5 * (x1 + x3);
        y13=_values->y(x13);
        x32=0.5 * (x3 + x2);
        y32=_values->y(x32);
        if(y3 > y13 && y3 > y32) {
          x1=x13;
          y1=y13;
          x2=x32;
          y2=y32;
        } else if(y13 > y32) {
          x2=x3;
          y2=y3;
          x3=x13;
          y3=y13;
        } else {
          x1=x13;
          y1=y13;
          x3=x32;
          y3=y32;
        }
      }
      if(isFound(y3, _y2)) {
        // two hidden roots are found, return the first one
        _x2=x3;
        _y2=y3;
        return true;
      }
      // maximum without root restore and go further
    }
    _x1=_x2;
    _y1=_y2;
    _x2=it.next(_x1);
  }
  // reaching the minimum slowness, tests if there's a last root
  if( ! it.atEnd(_x1, min, max)) {
    _x2=it.end(min, max);
    _y2=_values->y(_x2);
    if(isFound(_y1, _y2)) {
      return true;
    }
  }
  return false;
}

template <class RealType, class FunctionClass>
bool RootSolverTemplate<RealType, FunctionClass>::setInterval(RealType min, RealType max)
{
  _x1=min;
  _y1=_values->y(_x1);
  _x2=max;
  _y2=_values->y(_x2);
  return isFound(_y1, _y2);
}

template <class RealType, class FunctionClass>
void RootSolverTemplate<RealType, FunctionClass>::halving()
{
  // First halving
  RealType x3=0.5 * (_x1 + _x2);
  RealType y3=_values->y(x3);
  // Relative Convergence criteria
  RealType tmpAbs1;
  tmpAbs1=_x1 + _x2;
  if(tmpAbs1 < 0) tmpAbs1=-tmpAbs1;
  RealType maxDiff=tmpAbs1 * _precision;
  RealType minDiff=-maxDiff;
  for(register int i=1;i < ROOTSOLVER_MAX_ITERATIONS;i++ ) {
    if(y3==0.0) {
      // just move the X3 by one tenth of actuel range
      x3 -= (_x2 - _x1) * 0.1;
      // recalculate the Y3, it will always be different from 0
      y3=_values->y(x3);
    }
    if(isFound(_y1, y3)) {
      _x2=x3;
      _y2=y3;
    } else {
      _x1=x3;
      _y1=y3;
    }
    // Check for convergence with a relative criteria, if so exit
    RealType diff=_x1 - _x2;
    if(minDiff <= diff && diff <= maxDiff) {
      return;
    }
    x3=0.5 * (_x1 + _x2);
    y3=_values->y(x3);
  }
  App::stream() << "Halving maximum iteration reached" << endl;
}

template <class RealType, class FunctionClass>
void RootSolverTemplate<RealType, FunctionClass>::neville()
{
//#define neville_DEBUG
//#define neville_STAT
#if defined(neville_DEBUG) || defined(neville_STAT)
  printf("RootSolverTemplate::neville() @ precision=%lg\n", _precision);
#endif
  // Variables for polynomial fit of Neville
  int m=0;
  RealType p[ 20 ], nevY[ 20 ];
  // First halving
  RealType x3=0.5 * (_x1 + _x2);
  RealType y3=_values->y(x3);
  // Relative Convergence criteria
  RealType maxDiff=fabs(x3)*_precision;
  for (register int i=1; i<ROOTSOLVER_MAX_ITERATIONS; i++) {
    // define new bounds according to signs of y1 and y3
    // If y1 and y3 have different signs then the root is between
    // y1 and y3, thus forget (x2,y2) else forget (x1,y1).
    // this case is not implemented in Herrmann's code
    // In our case, it is very important to closely bracket the root
    // X1 and X2 must ALWAYS stay on the same sides of the root during all
    // Neville's iterations.
#ifdef neville_DEBUG
    //printf("1/_x1=%.20lf _y1=%.20lg 1/_x2=%.20lf _y2=%.20lg\n",1/_x1,_y1,1/_x2,_y2);
    printf("_x1=%.20lf _y1=%.20lg _x2=%.20lf _y2=%.20lg x3=%.20lf y3=%.20lg diff=%.20lg\n", _x1, _y1, _x2, _y2, x3, y3, fabs(_x1-_x2));
#endif
    if(y3==0.0) {
      // Move the x3 towards the furthest limit
      if(fabs(_x1-x3)>fabs(_x2-x3)) {
        x3=0.5*(_x1+x3);
      } else {
        x3=0.5*(_x2+x3);
      }
      // recalculate the Y3, it will always be different from 0
      y3=_values->y(x3);
#ifdef neville_DEBUG
      printf("y3 is null: x3=%.20lf y3=%.20lg\n", x3, y3);
#endif
    }
    if(isFound(_y1, y3)) {
      _x2=x3;
      _y2=y3;
    } else {
      _x1=x3;
      _y1=y3;
    }
    // Check for convergence with a relative criteria, if so exit
    if(fabs(_x1-_x2)<=maxDiff) {
#if defined(neville_DEBUG) || defined(neville_STAT)
      printf( "Number of iterations=%i \n", i);
      printf( "Neville poly m=%i \n", m);
      printf( "Convergence achieved\n" );
#endif
      return;
    }
    //  Check the slopes between x1,x3 and x3,x2
    //  If they are different, the function is strongly non linear,
    //  thus use only halving, not Neville
    if(isFound(_y1 - y3, y3 - _y2)) {
#if defined(neville_DEBUG) || defined(neville_STAT)
      printf( "Neville poly m=%i \n", m);
      printf( "Inverse slopes: " );
#endif
      x3=0.5* (_x1 + _x2);
      y3=_values->y(x3);
      m=0;
    } else {
      // Until 2009-02-03, this test was performed, found to be useless
      // If y1 and y2 differ by more than a factor of 100
      // use only halving to avoid poor behavior of polynomial fit

#if 0
      RealType ay1=_y1;
      if(ay1 < 0) ay1=-ay1;
      RealType ay2=_y2;
      if(ay2 < 0) ay2=-ay2;
      if(ay1 > 100 * ay2 || ay2 > 100 * ay1) {
#if defined(neville_DEBUG) || defined(neville_STAT)
        printf( "Neville poly m=%i \n", m);
        printf( "Factor of more than 100: " );
#endif
        x3=0.5* (x1 + x2);
        y3=_values->y(x3);
        m=0;
      } else {
#endif
#ifdef neville_DEBUG
      printf( "Neville poly m=%i \n", m);
#endif
      //    if(m>10) cout << "**** WARNING: m=" << m << endl;
      if(m>0) {
        p[ m + 1 ]=x3;
        nevY[ m + 1 ]=y3;
      } else {
        p[ 0 ]=_x1;
        nevY[ 0 ]=_y1;
        p[ 1 ]=_x2;
        nevY[ 1 ]=_y2;
      }
      // perform Neville iteration. Instead of generating y(x)
      // we interchange the x and y of formula to solve for x(y)
      // when y=0.
      // A new point is added to the pyramidal contruction.
      // The new point is always added at the base, adding a new
      // slice at each iteration, each time larger
      // As pyramide grows, the degree of the polynome increase
      // thus the number of point taken into account increase also.
      for(int j=m;j >= 0;j-- ) {
        RealType denom=nevY[ m + 1 ] - nevY[ j ];
        RealType tmpAbs1=nevY[ m + 1 ];
        if(tmpAbs1 < 0) tmpAbs1=-tmpAbs1;
        RealType tmpAbs2=denom;
        if(tmpAbs2 < 0) tmpAbs2=-tmpAbs2;
        if(tmpAbs2 < 1.0e-10 * tmpAbs1) {
#if defined(neville_DEBUG) || defined(neville_STAT)
          printf( "Neville poly m=%i \n", m);
          printf( "Denom too small (%lf), halving\n", denom);
#endif
          p[ 0 ]=0.5 * (_x1 + _x2);
          m=0;
          break;
        } else p[ j ]=(nevY[ m + 1 ] * p[ j ] - nevY[ j ] * p[ j + 1 ] )/denom;
      }
      x3=p[ 0 ];
      // Just move the X3 by one tenth of actual range towards the limit with the highest y
      RealType tmpAbs1=_y1<0 ? -_y1 : _y1;
      RealType tmpAbs2=_y2<0 ? -_y2 : _y2;
      if(tmpAbs1 > tmpAbs2)
        x3 -= 0.1 * (x3 - _x1);
      else
        x3 -= 0.1 * (x3 - _x2);
      //printf("Function estimation\n");
      m++;
      if(m > 18) m=1;
      // Make sure new estimate is inside the previous values
      // if not perform interval halving
      if(_x1 <= _x2) {
        if(x3 < _x1 || x3 > _x2) {
#if defined(neville_DEBUG) || defined(neville_STAT)
          printf( "Neville poly m=%i \n", m);
          printf( "New estimate out of range, halving\n" );
#endif
          x3=0.5* (_x1 + _x2);
          m=1;
        }
      } else {
        if(x3 < _x2 || x3 > _x1) {
#if defined(neville_DEBUG) || defined(neville_STAT)
          printf( "Neville poly m=%i \n", m);
          printf( "New estimate out of range, halving\n" );
#endif
          x3=0.5* (_x1 + _x2);
          m=1;
        }
      }
      y3=_values->y(x3);
    }
  }
  App::stream() << "Neville maximum iteration reached" << endl;
}

template <class RealType, class FunctionClass>
void RootSolverTemplate<RealType, FunctionClass>::newtonRaphson()
{
  // First halving
  RealType x3=0.5 * (_x1 + _x2);
  RealType y3=_values->y(x3);
  // Relative Convergence criteria
  RealType tmpAbs1;
  tmpAbs1=_x1 + _x2;
  if(tmpAbs1 < 0) tmpAbs1=-tmpAbs1;
  RealType maxDiff=tmpAbs1 * _precision;
  RealType minDiff=-maxDiff;
  for(register int i=1;i < ROOTSOLVER_MAX_ITERATIONS;i++ ) {
    if(y3==0.0) {
      // just move the X3 by one tenth of actuel range
      x3 -= (_x2 - _x1) * 0.1;
      // recalculate the Y3, it will always be different from 0
      y3=_values->y(x3);
    }
    if(isFound(_y1, y3)) {
      _x2=x3;
      _y2=y3;
    } else {
      _x1=x3;
      _y1=y3;
    }
    // Check for convergence with a relative criteria, if so exit
    RealType diff=_x1 - _x2;
    if(minDiff <= diff && diff <= maxDiff) {
      return;
    }
    // Compute the derivative estimation from the lowest y
    // Just move the X3 by one tenth of actual range towards the limit with the highest y
    RealType tmpAbs1=_y1<0 ? -_y1 : _y1;
    RealType tmpAbs2=_y2<0 ? -_y2 : _y2;
    if(tmpAbs1 > tmpAbs2) {
      x3=_x2 - _y2/_values->derivative(_x2);
      x3 -= 0.1 * (x3 - _x1);
    } else {
      x3=_x1 - _y1/_values->derivative(_x1);
      x3 -= 0.1 * (x3 - _x2);
    }
    // Make sure new estimate is inside the previous values
    // if not perform interval halving
    if(_x1 <= _x2) {
      if(x3 < _x1 || x3 > _x2) {
        x3=0.5* (_x1 + _x2);
      }
    } else {
      if(x3 < _x2 || x3 > _x1) {
        x3=0.5* (_x1 + _x2);
      }
    }
    y3=_values->y(x3);
  }
  App::stream() << "Newton-Raphson maximum iteration reached" << endl;
}

} // namespace QGpCoreTools

#endif // ROOTSOLVER_H