      real function eps(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,
*     eps(1.0) is known as the largest relative
*     spacing, beta**(1-t), and eps(-1.0) is the
*     smallest relative spacing, beta**(-t).
*
*     (14-Feb-1992)
      logical isInf, isNaN
      real x, beta, betain, temp

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

      beta = 2.0
      betain = 1.0/beta

      if (isInf(x)) then
*          x is +infinity or -infinity
           eps = x
      else if (isNaN(x)) then
*          x is a NaN
           eps = x
      else
           if (x .eq. 0.0) then
                eps = 1.0
           else
                eps = abs(x)
           end if
 10        temp = eps*betain
           temp = x + temp
           if (temp .ne. x) then
                eps = eps*betain
                go to 10
           end if
*          purify eps to eliminate all but high-order bit.
           if (isinf(x + eps)) then
                eps = x - eps
                eps = x - eps
           else
                eps = x + eps
                eps = eps - x
           end if
      end if
      end
