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

# include "zero_itp.h"

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

void zero_itp ( double f ( double ), double a, double b, double epsi, double k1, 
  double k2, int n0, bool verbose, double *z, double *fz, int *calls )

/******************************************************************************/
/*
  Purpose:

    zero_itp() seeks a zero of a function using the ITP algorithm.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    02 March 2024

  Author:

    John Burkardt

  Input:

    function f ( x ): the name of the user-supplied function.

    real a, b: the endpoints of the change of sign interval.

    real epsi: error tolerance between exact and computed roots.
    A reasonable value might be sqrt(eps).

    real k1: a parameter, with suggested value 0.2 / ( b - a ).

    real k2: a parameter, typically set to 2.

    int n0: a parameter that can be set to 0 for difficult problems,
    but is usually set to 1, to take more advantage of the secant method.

    bool verbose: if true, requests additional printed output.
*/
{
  double c;
  double delta;
  int nh;
  int nmax;
  double r;
  double s;
  double sigma;
  double xf;
  double xh;
  double xitp;
  double xt;
  double ya;
  double yb;
  double yitp;

  if ( b < a )
  {
    c = a;
    a = b;
    b = c;
  }

  ya = f ( a );
  yb = f ( b );
  if ( 0.0 < ya * yb )
  {
    printf ( "\n" );
    printf ( "zero_itp(): Fatal error!\n" );
    printf ( "  f(a) and f(b) have sign sign.\n" );
  }
/*
  Modify f(x) so that y(a) < 0, 0 < y(b);
*/
  if ( 0.0 < ya )
  {
    s = -1.0;
    ya = - ya;
    yb = - yb;
  }
  else
  {
    s = +1.0;
  }
  
  nh = ceil ( log2 ( ( b - a ) / 2.0 / epsi ) );
  nmax = nh + n0;

  *calls = 0;

  if ( verbose )
  {
    printf ( "\n" );
    printf ( "  User has requested additional verbose output.\n" );
    printf ( "  step   [a,    b]    x    f(x)\n" );
    printf ( "\n" );
  }

  while ( 2.0 * epsi < ( b - a ) )
  {
/*
  Calculate parameters.
*/
    xh = 0.5 * ( a + b );
    r = epsi * pow ( 2.0, nmax - *calls ) - 0.5 * ( b - a );
    delta = k1 * pow ( b - a, k2 );
/*
  Interpolation.
*/
    xf = ( yb * a - ya * b ) / ( yb - ya );
/*
  Truncation.
*/
    if ( 0 <= xh - xf )
    {
      sigma = +1;
    }
    else
    {
      sigma = -1;
    }

    if ( delta < fabs ( xh - xf ) )
    {
      xt = xf + sigma * delta;
    }
    else
    {
      xt = xh;
    }
/*
  Projection.
*/
    if ( fabs ( xt - xh ) <= r )
    {
      xitp = xt;
    }
    else
    {
      xitp = xh - sigma * r;
    }
/*
  Update the interval.
*/
    yitp = s * f ( xitp );

    if ( verbose )
    {
      printf ( "%d  [%g,%g]  f(%g)=%g\n", *calls, a, b, xitp, yitp );
    }

    if ( 0.0 < yitp )
    {
      b = xitp;
      yb = yitp;
    }
    else if ( yitp < 0.0 )
    {
      a = xitp;
      ya = yitp;
    }
    else
    {
      a = xitp;
      b = xitp;
      break;
    }
    *calls = *calls + 1;
  }

  *z = 0.5 * ( a + b );
  *fz = f ( *z );

  return;
}