SQRT(P**2 + Q**2)
PROBLEM
Compute SQRT(P**2+Q**2)
ALGORITHM
Moler and Morrison have described a fast algorithm for computing SQRT(P**2+Q**2). This algorithm has cubic convergence. It means that the result is accurate to 6.5 decimal digits after two iterations, to 20 digits after three iterations, and to 62 digits after four iterations.
IMPLEMENTATION
Unit:
internal function, external function without procedure statement
Parameters: numbers P, Q, a positive integer IterCount - number of iterations
Returns: value SQRT(P**2+Q**2) with IterCount**3 decimal digits
SQRTP2AQ2: procedure
parse arg P, Q, IterCount
numeric digits IterCount ** 3
P = ABS(P); Q = ABS(Q)
if P < Q
then do; W = P; P = Q; Q = W; end
if P = 0 then return Q
do IterCount
R = Q / P; R = R * R; R = R / (4 + R)
P = P + 2 * R * P; Q = Q * R
end
return P
|
Literature
Bentley J., Programming Pearls
CACM, December 1986, Vol. 29, No. 12, p. 1161