# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "zero_itp.hpp"

//****************************************************************************80

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 )

//****************************************************************************80
//
//  Purpose:
//
//    zero_itp() seeks a zero of a function using the ITP algorithm.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    03 March 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double f ( double x ): the user-supplied function.
//
//    double a, b: the endpoints of the change of sign interval.
//
//    double epsi: error tolerance between exact and computed roots.
//    A reasonable value might be sqrt(eps).
//
//    double k1: a parameter, with suggested value 0.2 / ( b - a ).
//
//    double 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.
//
//  Output:
//
//    double &z, &fz: the estimated root, and its function value.
//
//    int &calls: the number of function calls made.
//
{
  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 )
  {
    cout << "\n";
    cout << "zero_itp(): Fatal error!\n";
    cout << "  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 )
  {
    cout << "\n";
    cout << "  User has requested additional verbose output.\n";
    cout << "  step   [a,    b]    x    f(x)\n";
    cout << "\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 )
    {
      cout << calls 
           << "[" << a
           << "," << b
           << "] " << xitp
           << yitp << "\n";
    }

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