real function sqrt(x)
*     Cody-Waite implementation of sqrt(x)
*     [13-Nov-1990]

      integer and, e, intxp
      real f, setxp, x, y, z

      real onemme
      real Inf
*     real NaN

*     onemme = 1.0 - machine epsilon
      data onemme /z'3f7fffff'/
      data Inf    /z'7f800000'/
*     data NaN    /z'7fffffff'/

*     Generate a NaN for x <= 0.0, and Inf for x == Inf
*     and handle the special case of x == 0.0 so we do
*     not violate the assumptions that the arguments to
*     setxp() and intxp() are non-zero.

      if (x .eq. 0.0) then
           sqrt = 0.0
           return
      else if (x .lt. 0.0) then
           sqrt = 0.0
           sqrt = sqrt/sqrt
           return
      else if (x .ge. Inf) then
           sqrt = 0.0
           sqrt = 1.0/sqrt
           return
      else if (x .ne. x) then
           sqrt = 0.0
           sqrt = sqrt/sqrt
           return
      end if

      e = intxp(x)
      f = setxp(x,0)

*     y0 to 7.04 bits
      y = 0.41731 + 0.59016 * f

*     y1 to 15.08 bits
      z = y + f/y

*     y2 to 31.16 bits
      y = 0.25*z + f/z

      if (and(e,1) .ne. 0) then
         y = y * 0.70710 67811 86547 52440 08443 62104 E+00
         e = e + 1
      end if

*     Ensure 0.5 .le. y .lt. 1.0; otherwise,
*     setxp() will give wrong answer
      y = min(y,onemme)
      sqrt = setxp(y,e/2)

      end