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