*=======================================================================
* program main
* test program for root finding routine.
      program main
      implicit none
      double precision func0, func1, func2, func3, func4, func5, func6,
     &  a0, b0, xend
      integer nmbf
      external findroot, func0, func1, func2, func3, func4,
     &  func5, func6
      write (*,'(2x,a)') 'Bisection method'
      write (*,'(2x,a8,2x,a16,2x,a16,2x,a6)')
     &  'function', 'xend', 'error', 'nmbf'
      a0 = 1
      b0 = 2
      call findroot_bisect (func0, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func0', xend, xend-sqrt(2.0d0), nmbf
      a0 = 1
      b0 = 2
      call findroot_bisect (func1, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func1', xend, xend-sqrt(2.0d0), nmbf
      a0 = 1
      b0 = 2
      call findroot_bisect (func2, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func2', xend, xend-sqrt(2.0d0), nmbf
      a0 = 1
      b0 = 2
      call findroot_bisect (func3, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func3', xend, xend-sqrt(2.0d0), nmbf
      a0 = 1
      b0 = 2
      call findroot_bisect (func4, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func4', xend, xend-sqrt(2.0d0), nmbf
      a0 = 1
      b0 = 2
      call findroot_bisect (func5, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func5', xend, xend-sqrt(2.0d0), nmbf
      a0 = -1
      b0 = 2
      call findroot_bisect (func6, a0, b0, xend, nmbf)
      write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
     &  'func6', xend, xend-0.0d0, nmbf
      stop
      end
*=======================================================================
* subroutine findroot_bisect
      subroutine findroot_bisect (f, a0, b0, xend, nmbf)
      implicit none
      double precision f, a0, b0, xend
      integer nmbf
      external f
*     ------------------------------------------------------------------
      double precision signum, bisect, a, b, c, fa, fb, fc
      external signum, bisect
      a = a0
      b = b0
      fa = f(a)
      fb = f(b)
      nmbf = 2
      if (fa.eq.0) then
       c = a
       fc = fa
      else if (fb.eq.0) then
       c = b
       fc = fb
      else
       if (signum(fa).eq.signum(fb)) then
        stop 'error in findroot_bisect -- no bracket for root'
       endif
       c = bisect(a,b)
       fc = f(c)
       nmbf = nmbf+1
      endif
      do while (c.ne.a.and.c.ne.b)
       if (fc.eq.0) then
        a = c
        fa = fc
        b = c
        fb = fc
       else if (signum(fc).eq.signum(fa)) then
        a = c
        fa = fc
        c = bisect(a,b)
        fc = f(c)
        nmbf = nmbf+1
       else if (signum(fc).eq.signum(fb)) then
        b = c
        fb = fc
        c = bisect(a,b)
        fc = f(c)
        nmbf = nmbf+1
       else
        stop 'error in findroot_bisect -- f generates NaN?'
       endif
      enddo
      xend = c
      return
      end
*=======================================================================
      double precision function bisect(a,b)
      implicit none
      double precision a, b
*     ------------------------------------------------------------------
*     Modified bisection routine.  If a and b are both non-zero, of the
*     same sign, and not too far apart on a logarithmic scale, returns
*     c=(a+b)/2.  Other cases ... see the code.  The intent is that
*     bisect(a,b) divides the interval so that there are about equal
*     number of floating point values in either half.  Some preference
*     towards zero.
      double precision signum, eps, tiny, t0, t1
      external signum
*     let eps be IEEE machine precision
      eps = 2.0d0**(-52)
*     set tiny to IEEE smallest positive normal number
      tiny = 2.0d0**(-1022)
      if (signum(a)*signum(b).lt.0) then
       bisect = 0
      else if (signum(a)*signum(b).gt.0) then
       t0 = min(abs(a),abs(b))
       t1 = max(abs(a),abs(b))
       if (t1.lt.4*t0.or.t1.lt.4*tiny) then
        bisect = (a+b)/2
       else
        bisect = signum(a)*sqrt(t0)*sqrt(t1)
       endif
      else if (a.eq.0.or.b.eq.0) then
       t0 = a+b
       if (abs(t0).ge.1/sqrt(eps)) then
        bisect = signum(t0)
       else if (abs(t0).ge.sqrt(eps)) then
        bisect = t0*sqrt(eps)
       else if (abs(t0).ge.sqrt(tiny)) then
        bisect = t0*abs(t0)
       else if (abs(t0).ge.4*tiny) then
        bisect = signum(t0)*tiny
       else if (abs(t0).gt.eps*tiny) then
        bisect = signum(t0)*eps*tiny
       else
        bisect = 0
       endif
      else
*      must be NaN
       bisect = (a+b)/2
      endif
      return
      end
*=======================================================================
* double precision function signum(x)
* returns the sign of x; returns 0 if x.eq.0.
      double precision function signum(x)
      implicit none
      double precision x
      if (x.lt.0) then
       signum = -1
      else if (x.eq.0) then
       signum = 0
      else if (x.gt.0) then
       signum = 1
      else
*      (must be that x is Not-a-Number)
       signum = x
      endif
      return
      end
*=======================================================================
* double precision function func0(x)
* test function for rootfinding routine.
* func0 is smooth and has a simple root at x=sqrt(2.0d0).
      double precision function func0(x)
      implicit none
      double precision x
      func0 = (x-sqrt(2.0d0))*sqrt(1+x**2)
      return
      end
*=======================================================================
* double precision function func1(x)
* test function for rootfinding routine.
* func1 is smooth and has a triple root at x=sqrt(2.0d0).
      double precision function func1(x)
      implicit none
      double precision x
      func1 = (x-sqrt(2.0d0))**3
      return
      end
*=======================================================================
* double precision function func2(x)
* test function for rootfinding routine.
* func2 is smooth and has a simple root at x=sqrt(2.0d0); the
* asymptotic behavior may pose problems for a Newton-like method.
      double precision function func2(x)
      implicit none
      double precision x
      func2 = atan(x-sqrt(2.0d0))
      return
      end
*=======================================================================
* double precision function func3(x)
* test function for rootfinding routine.
* func3 is smooth and has a simple root at x=sqrt(2.0d0); the
* asymptotic behavior may pose problems for a Newton-like method.
      double precision function func3(x)
      implicit none
      double precision x
      func3 = atan(256*(x-sqrt(2.0d0)))
      return
      end
*=======================================================================
* double precision function func4(x)
* test function for rootfinding routine.
* func4 is discontinuous with sign change at x=sqrt(2.0d0).
      double precision function func4(x)
      implicit none
      double precision x
      if (x.lt.0.or.x*x.lt.2) then
       func4 = -1
      else
       func4 = 1
      endif
      return
      end
*=======================================================================
* double precision function func5(x)
* test function for rootfinding routine.
* func5 has a vertical asymptote and sign change at x=sqrt(2.0d0).
      double precision function func5(x)
      implicit none
      double precision x
      double precision t
      t = (x-sqrt(2.0d0))*sqrt(1+x**2)
      if (t.ne.0) then
       func5 = 1/t
      else
       func5 = 0
      endif
      return
      end
*=======================================================================
* double precision function func6(x)
* test function for rootfinding routine.
* func6 is discontinous with sign change at 0.
      double precision function func6(x)
      implicit none
      double precision x
      if (x.lt.0) then
       func6 = -1
      else
       func6 = 1
      endif
      return
      end