using Printf

function zero_itp( f, a, b, epsi, k1, k2, n0, verbose )

#*****************************************************************************80
#
## zero_itp() seeks the root of a function using Interpolate/Truncate/Project.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 March 2024
#
#  Author:
#
#    Original Julia code by TheLateKronos.
#    This version by 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.
#
#  Output:
#
#    real z, fz: the estimated root and its function value.
#
#    integer ncalls: the number of function calls.
#

#
#  Force a < b.
#
  if ( b < a )
    c = a
    a = b
    b = c
  end
#
#  Require f(a) and f(b) of opposite sign.
#
  ya = f( a )
  yb = f( b )

  if ( 0 < ya * yb )
    print( "" )
    print( "zero_itp(): Fatal error!" )
    print( "  f(a) and f(b) have same sign." )
  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

  n1_2 = ceil( Int, log2( ( b - a ) / 2.0 / epsi ) )
  nmax = n1_2 + n0

  calls = 0

  if ( verbose )
    print( "\n" )
    print( "  User has requested additional verbose output.\n" )
    print( "  step   [a,    b]    x    f(x)\n" )
    print( "\n" )
  end

  while ( 2.0 * epsi < ( b - a ) )
#
#  Calculate parameters.
#
    xh = ( a + b ) / 2.0
    r = epsi * 2.0 ^ ( nmax - calls ) - ( b - a ) / 2.0
    delta = k1 * ( b - a ) ^ k2
#
#  Interpolation:
#
    xf = ( yb * a - ya * b ) / ( yb - ya )
#
#  Truncate
#
    sigma = sign( xh - xf )

    if ( delta < abs( xh - xf ) )
      xt = xf + sigma * delta
    else
      xt = xh
    end
#
#  Project
#
    if ( abs( xt - xh ) <= r )
      xitp = xt
    else
      xitp = xh - sigma * r
    end
#
#  Update Interval.
#
    yitp = s * f( xitp )

    if ( verbose )
      @printf( "%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
    end

    calls = calls + 1

  end

  z = ( a + b ) / 2.0
  fz = f( z )

  return z, fz, calls
end