function [ z, fz, calls ] = zero_itp ( f, a, b, epsi, k1, k2, n0, verbose )

%*****************************************************************************80
%
%% 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.
%
%    real k1: a parameter, with suggested value 0.2 / ( b - a ).
%
%    real k2: a parameter, typically set to 2.
%
%    integer 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.
%
%    logical verbose: if true, requests additional printed output.
%
  if ( b < a )
    c = a;
    a = b;
    b = c;
  end

  ya = f ( a );
  yb = f ( b );
  if ( 0.0 < ya * yb )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'zero_itp(): Fatal error!\n' );
    fprintf ( 1, '  f(a) and f(b) have sign sign.\n' );
    error ( 'zero_itp(): Fatal error!' );
  end
%
%  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;
  end

  if ( nargin < 7 )
    n0 = 1;
  end

  if ( nargin < 6 )
    k2 = 2.0;
  end

  if ( nargin < 5 )
    k1 = 1.0 / 5.0 / ( b - a );
  end

  if ( nargin < 4 )
    epsi = sqrt ( eps );
  end
  
  nh = ceil ( log2 ( ( b - a ) / 2.0 / epsi ) );
  nmax = nh + n0;

  calls = 0;

  if ( verbose )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User has requested additional verbose output.\n' );
    fprintf ( 1, '  step   [a,    b]    x    f(x)\n' );
    fprintf ( 1, '\n' );
  end

  while ( 2.0 * epsi < ( b - a ) )
%
%  Calculate parameters.
%
    xh = 0.5 * ( a + b );
    r = epsi * 2.0 ^ ( nmax - calls ) - 0.5 * ( b - a );
    delta = k1 * ( b - a ) ^ k2;
%
%  Interpolation.
%
    xf = ( yb * a - ya * b ) / ( yb - ya );
%
%  Truncation.
%
    sigma = sign ( xh - xf );
    if ( delta < abs ( xh - xf ) )
      xt = xf + sigma * delta;
    else
      xt = xh;
    end
%
%  Projection.
%
    if ( abs ( xt - xh ) <= r )
      xitp = xt;
    else
      xitp = xh - sigma * r;
    end
%
%  Update the interval.
%
    yitp = s * f ( xitp );

    if ( verbose )
      fprintf ( 1, '%d  [%g,%g]  f(%g)=%g\n', calls, a, b, xitp, yitp );
    end

    if ( 0.0 < yitp )
      b = xitp;
      yb = yitp;
    elseif ( yitp < 0.0 )
      a = xitp;
      ya = yitp;
    else
      a = xitp;
      b = xitp;
      break;
    end

    calls = calls + 1;

  end

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

  return
end

