real function  sqrt(xs)
*     Cody-Waite implementation of sqrt(x)
*     Intermediate computation all in double precision
*     [14-Nov-1990]

      integer and, e, dintxp
      real xs
      double precision dsetxp, f, x, y, z

      double precision onemme
      double precision Inf
*     double precision NaN

*     onemme = 1.0 - machine epsilon
      data onemme /z'3fefffffffffffff'/
      data Inf    /z'7ff0000000000000'/
*     data NaN    /z'7fffffffe0000000'/

      x = dble(xs)

*     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.0D+00) then
           sqrt = 0.0D+00
           return
      else if (x .lt. 0.0D+00) then
           sqrt = 0.0D+00
           sqrt = sqrt/sqrt
           return
      else if (x .ge. Inf) then
           sqrt = 0.0D+00
           sqrt = 1.0D+00/sqrt
           return
      else if (x .ne. x) then
           sqrt = 0.0D+00
           sqrt = sqrt/sqrt
           return
      end if

      e = dintxp(x)
      f = dsetxp(x,0)

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

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

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

*     y3 to 63.32 bits
      y = 0.5D+00*(y + f/y)

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

*     Ensure 0.5 .le. y .lt. 1.0; otherwise,
*     dsetxp() will give wrong answer
      y = min(max(0.5D+00,y),onemme)

*     Finally, reduce double precision result to single precision
      sqrt = sngl(dsetxp(y,e/2))

      end