# include <cmath>
# include <cstdlib>

using namespace std;

# include "zero_chandrupatla.hpp"

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

void zero_chandrupatla ( double f ( double x ), double x1, double x2, 
  double &xm, double &fm, int &calls )

/******************************************************************************/
//
//  Purpose:
//
//    zero_chandrupatla() seeks a zero of a function using Chandrupatla's algorithm.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 March 2024
//
//  Author:
//
//    Original QBASIC version by Tirupathi Chandrupatla.
//    This version by John Burkardt.
//
//  Reference:
//
//    Tirupathi Chandrupatla,
//    A new hybrid quadratic/bisection algorithm for finding the zero of a
//    nonlinear function without using derivatives,
//    Advances in Engineering Software,
//    Volume 28, Number 3, pages 145-149, 1997.
//
//  Input:
//
//    double f ( double x ): the name of the user-supplied function.
//
//    double a, b: the endpoints of the change of sign interval.
//
//  Output:
//
//    double &z, &fz: the estimated root and its function value.
//
//    int &calls: the number of function calls.
//
{
  double al;
  double a;
  double b;
  double c;
  double d;
  const double delta = 0.00001;
  const double epsilon = 1.0E-10;
  double f0;
  double f1;
  double f2;
  double f3;
  double fh;
  double fl;
  double ph;
  double t;
  double tl;
  double tol;
  double x0;
  double x3;
  double xi;

  f1 = f ( x1 );
  f2 = f ( x2 );
  calls = 2;

  t = 0.5;

  while ( true )
  {
    x0 = x1 + t * ( x2 - x1 );
    f0 = f ( x0 );
    calls = calls + 1;
//
//  Arrange 2-1-3: 2-1 Interval, 1 Middle, 3 Discarded point.
//
    if ( ( 0 < f0 ) ==  ( 0 < f1 ) )
    {
      x3 = x1;
      f3 = f1;
      x1 = x0;
      f1 = f0;
    }
    else
    {
      x3 = x2;
      f3 = f2;
      x2 = x1;
      f2 = f1;
      x1 = x0;
      f1 = f0;
    }
//
//  Identify the one that approximates zero.
//
    if ( fabs ( f2 ) < fabs ( f1 ) )
    {
      xm = x2;
      fm = f2;
    }
    else
    {
      xm = x1;
      fm = f1;
    }

    tol = 2.0 * epsilon * fabs ( xm ) + 0.5 * delta;
    tl = tol / fabs ( x2 - x1 );

    if ( 0.5 < tl || fm == 0.0 )
    {
      break;
    }
//
//  If inverse quadratic interpolation holds, use it.
//
    xi = ( x1 - x2 ) / ( x3 - x2 );
    ph = ( f1 - f2 ) / ( f3 - f2 );
    fl = 1.0 - sqrt ( 1.0 - xi );
    fh = sqrt ( xi );

    if ( fl < ph && ph < fh )
    {
      al = ( x3 - x1 ) / ( x2 - x1 );
      a = f1 / ( f2 - f1 );
      b = f3 / ( f2 - f3 );
      c = f1 / ( f3 - f1 );
      d = f2 / ( f3 - f2 );
      t = a * b + c * d * al;
    }
    else
    {
      t = 0.5;
    }
//
//  Adjust T away from the interval boundary.
//
    t = fmax ( t, tl );
    t = fmin ( t, 1.0 - tl );
  }

  return;
}