// RootFindLinear.cpp

#include "RootFindLinear.h"
#include <cmath>
#include <cassert>

#ifndef DEFECT_CMATH_NOT_STD
using std::fabs;
#endif

//**********************************************************************

// constructor

RootFindLinear::RootFindLinear() : _max_iter(100) { }

//**********************************************************************

// destructor

RootFindLinear::~RootFindLinear() { }

//**********************************************************************

// Find the root.

StatusDouble RootFindLinear::solve(double x1, double x2) const {

  // Value to return for failure.
  double xfail = 0.0;

  // Save limits.
  // We require solution to be in the range (xmin,xmax).
  double xmin = x1;
  double xmax = x2;

  if ( xmin >= xmax ) return StatusDouble(INVALID_LIMITS,xfail);

  // Evaluate the starting difference.
  double adif = fabs(x2-x1);

  // Set the precision.
  // Solution is said to be close if it changes by less than this value.
  double xclose = 1.e-14 * fabs(x2-x1);
  double yclose = 1.e-14;

  StatusDouble sd1 = evaluate(x1);
  if ( sd1.status() ) return StatusDouble(INVALID_FUNCTION_CALL,xfail);
  double f1 = sd1.value();
  StatusDouble sd2 = evaluate(x2);
  if ( sd2.status() ) return StatusDouble(INVALID_FUNCTION_CALL,xfail);
  double f2 = sd2.value();

  // Loop.
  for ( int count=0; count<_max_iter; ++count ) {
    // interpolate for new solution
    assert( f1 != f2 );
    if ( f1 == f2 ) return StatusDouble(OUT_OF_RANGE,xfail);
    double x0 = (f2*x1 - f1*x2) / (f2 - f1);
    if ( x0 < xmin ) return StatusDouble(OUT_OF_RANGE,xfail);
    if ( x0 > xmax ) return StatusDouble(OUT_OF_RANGE,xfail);
    // Evaluate the function at the new solution.
    StatusDouble sd0 = evaluate(x0);
    if ( sd0.status() ) return StatusDouble(INVALID_FUNCTION_CALL,xfail);
    double f0 = sd0.value();
    // Replace the most distant guess with x0.
    double adif1 = fabs(x1-x0);
    double adif2 = fabs(x2-x0);
    if ( adif1 < adif2 ) {
      adif = adif1;
      x2 = x0;
      f2 = f0;
    }
    else {
      adif = adif2;
      x1 = x0;
      f1 = f0;
    }
    // Return if solution is found.
    if ( f0 == 0.0 ) return StatusDouble(0,x0);
    // If we are close in y, exit.
    if ( fabs(f2-f1) < yclose ) return StatusDouble(0,x0);
    // If we are close in x, exit.
    if ( adif < xclose ) return StatusDouble(0,x0);
  }

  // Too many iterations.
  return StatusDouble(TOO_MANY_ITERATIONS,xfail);

} 

//**********************************************************************