double precision function deps(x)
*     (Machine Epsilon)
*
*     Return the machine epsilon, the smallest
*     positive number such that (x + epsilon) .ne.
*     x.  This code works for x = 0.0, NaN, and
*     Inf as well.
*
*     Note that by supplying negative x values,
*     you can determine the smallest positive
*     number such that (x - epsilon) .ne. x.
*
*     For base beta and t-digit significands,
*     deps(1.0) is known as the largest relative
*     spacing, beta**(1-t), and deps(-1.0) is the
*     smallest relative spacing, beta**(-t).
*
*     (14-Feb-1992)
      logical disInf, disNaN
      double precision x, beta, betain, temp

*     WARNING: beta is machine dependent!
*     parameter (beta = 2.0)
*     parameter (betain = 1.0/beta)

      beta = 2.0d+00
      betain = 1.0d+00/beta

      if (disInf(x)) then
*          x is +infinity or -infinity
           deps = x
      else if (disNaN(x)) then
*          x is a NaN
           deps = x
      else
           if (x .eq. 0.0d+00) then
                deps = 1.0d+00
           else
                deps = dabs(x)
           end if
 10        temp = deps*betain
           temp = x + temp
           if (temp .ne. x) then
                deps = deps*betain
                go to 10
           end if
*          purify deps to eliminate all but high-order bit.
           if (disInf(x + deps)) then
                deps = x - deps
                deps = x - deps
           else
                deps = x + deps
                deps = deps - x
           end if
      end if
      end