double precision function  dsqrt(x)
*     Cody-Waite implementation of dsqrt(x)
*     [14-Nov-1990]

      integer and, e, dintxp
      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'/

*     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
           dsqrt = 0.0D+00
           return
      else if (x .lt. 0.0D+00) then
           dsqrt = 0.0D+00
           dsqrt = dsqrt/dsqrt
           return
      else if (x .ge. Inf) then
           dsqrt = 0.0D+00
           dsqrt = 1.0D+00/dsqrt
           return
      else if (x .ne. x) then
           dsqrt = 0.0D+00
           dsqrt = dsqrt/dsqrt
           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(y,onemme)
      dsqrt = dsetxp(y,e/2)

      end