# include <math.h>
# include <stdbool.h>
# include <stdlib.h>

# include "zero_chandrupatla.h"

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

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:

    18 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;
}