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