itp Subroutine

    subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)

    implicit none

    class(itp_solver),intent(inout) :: me
    real(wp),intent(in)    :: ax      !! left endpoint of initial interval
    real(wp),intent(in)    :: bx      !! right endpoint of initial interval
    real(wp),intent(in)    :: fax     !! `f(ax)`
    real(wp),intent(in)    :: fbx     !! `f(ax)`
    real(wp),intent(out)   :: xzero   !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
    real(wp),intent(out)   :: fzero   !! value of `f` at the root (`f(xzero)`)
    integer,intent(out)    :: iflag   !! status flag (`0`=root found, `-2`=max iterations reached)

    real(wp) :: a,b,ya,yb,x12,r,d,xf,bma,sigma,xt,xitp,yitp,term,aprev,bprev,denom
    integer :: n12, nmax
    integer :: j !! iteration counter
    logical :: root_found !! convergence in x
    logical :: fail !! if we can't interpolate
    real(wp) :: factor !! factor to ensure f(ax)<f(bx)

    real(wp),parameter :: log2 = log(2.0_wp)

    ! initialize:
    if (fax < fbx) then  ! a requirement for the algorithm ?
        factor = 1.0_wp
    else
        ! flip the sign of the function to ensure fax < fbx
        ! note that a < b is already enforced by the calling routine.
        factor = -1.0_wp
    end if
    a  = ax
    b  = bx
    ya = factor*fax
    yb = factor*fbx
    iflag = 0
    term = (b-a)/(2.0_wp*me%rtol)
    n12 = ceiling ( log(abs(term)) / log2 ) ! ceiling(log2(term))
    nmax = n12 + me%n0
    aprev = huge(1.0_wp) ! initialize to unusual values
    bprev = huge(1.0_wp) !

    ! main loop
    do j = 0, min(nmax,me%maxiter)

        x12 = bisect(a,b)
        bma = b - a
        ! note: protect for r<0 as mentioned in paper
        r = max(me%rtol * 2.0_wp ** (nmax-j) - bma/2.0_wp, 0.0_wp)
        d = me%k1 * bma**me%k2

        ! interpolation:
        denom = yb - ya
        fail = abs(denom) <= tiny(1.0_wp)  ! check for divide by zero

        if (.not. fail) then
            xf = (yb*a - ya*b) / denom

            ! truncation:
            sigma = sign(1.0_wp, x12 - xf)
            if (d <= abs(x12-xf)) then
                xt = xf + sigma * d
            else
                xt = x12
            end if

            ! projection:
            if (abs(xt - x12) <= r) then
                xitp = xt
            else
                xitp = x12 - sigma * r
            end if

            ! updating interval:
            yitp = factor*me%f(xitp)
            if (me%solution(xitp,yitp,xzero,fzero)) return
            if (yitp > 0.0_wp) then
                b = xitp
                yb = yitp
            elseif (yitp < 0.0_wp) then
                a = xitp
                ya = yitp
            else
                a = xitp
                b = xitp
            end if

        end if

        if (fail .or. (a==aprev .and. b==bprev)) then
            ! if the interval hasn't changed, then it is stuck.
            ! [this can happen in the test cases].
            ! So just do a bisection
            xitp = bisect(a,b)
            yitp = factor*me%f(xitp)
            if (me%solution(xitp,yitp,xzero,fzero)) return
            if (ya*yitp<0.0_wp) then
                ! root lies between a and xitp
                b = xitp
                yb = yitp
            else
                ! root lies between xitp and b
                a = xitp
                ya = yitp
            end if
        end if
        aprev = a
        bprev = b

        ! check for convergence:
        root_found = me%converged(a,b)
        if (root_found .or. j==me%maxiter) then
            call choose_best(a,b,ya,yb,xzero,fzero)
            if (.not. root_found) iflag = -2  ! max iterations reached
            exit
        end if

    end do

    end subroutine itp