real function sqrt(x) ************************************************************************ * (Square Root) * This function computes the square root of its argument via * a cubically-convergent iterative algorithm published by * * Cleve Moler and Donald Morrison, "Replacing Square Roots by * Pythagorean Sums", IBM J. Research and Development, 27, * 577-581 (1983). * * Augustin A. Dubrulle, "A Class of Numerical Methods for the * Computation of Pythagorean Sums", IBM J. Research and * Development, 27, 582-589 (1983). * * The algorithm has been augmented by checks for negative, NaN, * and Infinity arguments. These exceptional conditions generate * run-time NaN (negative or NaN x) or Infinity values that can be * trapped if necessary. * * This algorithm is ABSYMAL for |x| << 1; for x = 0, after 10,000 * iterations, it had converged only to p = 1.e-6. * (12-Feb-1994) ************************************************************************ real ZERO parameter (ZERO = 0.0e+00) real ONE parameter (ONE = 1.0e+00) real TWO parameter (TWO = 2.0e+00) real FOUR parameter (FOUR = 4.0e+00) real p, r, s, x logical isNaN, isInf integer n if (x .lt. ZERO) then sqrt = ZERO sqrt = sqrt / sqrt else if (isNaN(x)) then sqrt = x + x else if (isInf(x)) then sqrt = x + x else n = 0 p = ONE r = x - ONE 10 if ((ONE + r) .ne. ONE) then s = r / (FOUR + r) p = p + TWO * s * p r = r * (s/(ONE + TWO * s))**2 n = n + 1 if (n .gt. 5) print *,'x =',x,' s = ',s,' p =',p, x ' r =',r if (n .gt. 10) then print *, x '---------- mmsqrt() not converging ----------' else go to 10 endif endif sqrt = p endif end