*=======================================================================
* 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,a8,2x,a16,2x,a16,2x,a6)')
     &  'function', 'xend', 'error', 'nmbf'
      a0 = 1
      b0 = 2
      call findroot (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 (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 (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 (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 (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 (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 (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
* Inputs:
*   f: double precision function of one double precision argument.
*   a0, b0: double precision bounding bracket for the root.
*   Must have signum(f(a0))*signum(f(b0)).le.0.
* Returns: 
*   xend: root of f in the interval defined by (a0,b0).
*   nmbf: number of calls to f made for this call to findroot.
      subroutine findroot (f, a0, b0, xend, nmbf)
      implicit none
      double precision f, a0, b0, xend
      integer nmbf
      external f
*     ------------------------------------------------------------------
      ...
      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