subroutine zero_itp ( f, a, b, epsi, k1, k2, n0, verbose, z, 
     &  fz, calls )

c*********************************************************************72
c
cc zero_itp() seeks a zero of a function using the ITP algorithm.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 March 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    function f ( x ): the name of the user-supplied function.
c
c    double precision a, b: the endpoints of the change of sign interval.
c
c    double precision epsi: error tolerance between exact and computed roots.
c
c    double precision k1: a parameter, with suggested value 0.2 / ( b - a ).
c
c    double precision k2: a parameter, typically set to 2.
c
c    integer n0: a parameter that can be set to 0 for difficult problems,
c    but is usually set to 1, to take more advantage of the secant method.
c
c    logical verbose: if true, requests additional printed output.
c
c  Output:
c
c    double precision z, fz: the estimated root and its function value.
c
c    integer ncalls: the number of function calls.
c
      implicit none

      double precision a
      double precision b
      double precision c
      integer calls
      double precision delta
      double precision epsi
      double precision f
      external f
      double precision fz
      double precision k1
      double precision k2
      integer n0
      integer nh
      integer nmax
      double precision r
      double precision r8_log_2
      double precision s
      double precision sigma
      logical verbose
      double precision xf
      double precision xh
      double precision xitp
      double precision xt
      double precision ya
      double precision yb
      double precision yitp
      double precision z

      if ( b < a ) then
        c = a
        a = b
        b = c
      end if

      ya = f ( a )
      yb = f ( b )
      if ( 0.0D+00 < ya * yb ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'zero_itp(): Fatal errorc'
        write ( *, '(a)' ) '  f(a) and f(b) have sign sign.'
        stop ( 1 )
      end if
c
c  Modify f(x) so that y(a) < 0, 0 < y(b)
c
      if ( 0.0D+00 < ya ) then
        s = -1.0D+00
        ya = - ya
        yb = - yb
      else
        s = +1.0D+00
      end if
  
      nh = ceiling ( r8_log_2 ( ( b - a ) / 2.0D+00 / epsi ) )
      nmax = nh + n0

      calls = 0

      if ( verbose ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) '  User requested additional verbose output.'
        write ( *, '(a)' ) '  step   [a,    b]    x    f(x)'
        write ( *, '(a)' ) ''
      end if

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

        if ( verbose ) then
          write ( *, '(2x,i4,2x,f10.6,2x,f10.6,2x,f10.6,2x,g14.6)' ) 
     &      calls, a, b, xitp, yitp
        end if

        if ( 0.0D+00 < yitp ) then
          b = xitp
          yb = yitp
        else if ( yitp < 0.0D+00 ) then
          a = xitp
          ya = yitp
        else
          a = xitp
          b = xitp
          exit
        end if

        calls = calls + 1

      end do

      z = 0.5D+00 * ( a + b )
      fz = f ( z )

      return
      end
      function r8_log_2 ( x )

c*********************************************************************72
c
cc r8_log_2() returns the logarithm base 2 of an R8.
c
c  Discussion:
c
c    value = Log ( |X| ) / Log ( 2.0 )
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 August 2002
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision X, the number whose base 2 logarithm is desired.
c    X should not be 0.
c
c  Output:
c
c    double precision R8_LOG_2, the logarithm base 2 of the absolute
c    value of X.  It should be true that |X| = 2^R8_LOG_2.
c
      implicit none

      double precision r8_log_2
      double precision x

      if ( x == 0.0D+00 ) then
        r8_log_2 = - huge ( x )
      else
        r8_log_2 = log ( abs ( x ) ) / log ( 2.0D+00 )
      end if

      return
      end