-- © Copyright 1999 by John Halleck-- All rights reservedpackagebodyGeneric_Pythagorean_Sumsis-- Compute the Pythagorean Sum (I.E. what is the distance given Delta X, and-- DeltaY)-- Traditionally, this is done as sqrt(x*x+y*y), but this has the problem that-- the squared quantity can overflow long before the result would have.-- Cleve Moler and Donald Morrison discovered (as documented in their article-- in the IBM Journal of Research and Development, VOL27 No 6 NOV 83) that-- this could actually computed more directly.-- Advantages of this algorithm:-- It never squares X or Y, so it never overflows unless the-- answer does.-- It is very numericly stable.-- It is very very very fast.-- This code started life in MASM Assembley Language on a Univac.-- That code was then ported to Fortran IV, on a Univac.-- That code was then ported to C.-- That code was then used as the base of this port.-- SO.. it has a lot of inherited uglies.functionP_Sum (X, Y : Float_Type)returnFloat_TypeisGuess, Error : Float_Type; Temp : Float_Type; Delta_X, Delta_Y : Float_Type; R, S : Float_Type;begin-- First approximation.Delta_X :=abs(X); Delta_Y :=abs(Y); Guess := Float_Type'Max (Delta_X, Delta_Y); Error := Float_Type'Min (Delta_X, Delta_Y);-- The following code is the loop unrolling of:-- Temp := Error / Guess;-- R := Temp * Temp;-- S := R / (R + 4.0);-- Guess := Guess + 2.0 * S * Guess;-- Error := Error * S;-- The loop is unrolled because checking for convergence-- takes more time than the computation.-- 1 iteration gives 1 significant digits in the result.-- 2 iterations give 5 significant digits,-- 3 iterations give 20 significant digits,-- 4 iterations give 62 significant digits.-- Etc. (The number of significant decimal digits roughly triples-- with each iteration.)-- The loop invariant is that SQRT(GUESS**2+ERROR**2) is the Length-- with each iteration dramaticly shrinking ERROR.-- After this section we have one significant figure in the resultTemp := Error / Guess; R := Temp * Temp; S := R / (R + 4.0); Guess := Guess + 2.0 * S * Guess;-- After this next section we have 6 significant digits in the resultError := Error * S; Temp := Error / Guess; R := Temp * Temp; S := R / (R + 4.0); Guess := Guess + 2.0 * S * Guess;ifFloat_Type'Digits < 5thenreturnGuess;endif;-- After this section we have 20 significant digits in the result.Error := Error * S; Temp := Error / Guess; R := Temp * Temp; S := R / (R + 4.0); Guess := Guess + 2.0 * S * Guess;ifFloat_Type'Digits < 20thenreturnGuess;endif;-- After this section we have 62 significant digits in the result.Error := Error * S; Temp := Error / Guess; R := Temp * Temp; S := R / (R + 4.0); Guess := Guess + 2.0 * S * Guess;ifFloat_Type'Digits < 62thenreturnGuess;endif;-- After this section we have at least 120 significant digits.Error := Error * S; Temp := Error / Guess; R := Temp * Temp; S := R / (R + 4.0); Guess := Guess + 2.0 * S * Guess;returnGuess;endP_Sum;endGeneric_Pythagorean_Sums;

This page is http://www.cc.utah.edu/~nahaj/ada/math/generic_pythagorean_sums.adb.html

© Copyright 2000 by John Halleck, All Rights Reserved.

This snapshot was last modified on November 28th, 2000