# Generic_pythagorean_sums (generic_pythagorean_sums.adb): (Specification)

```--  © Copyright 1999 by John Halleck

package body Generic_Pythagorean_Sums is

--  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
--    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.

function P_Sum (X, Y : Float_Type) return Float_Type is

Guess, 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 result
Temp  := 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 result
Error := Error * S;
Temp  := Error / Guess;
R     := Temp * Temp;
S     := R / (R + 4.0);
Guess := Guess + 2.0 * S * Guess;

if Float_Type'Digits < 5 then return Guess; end if;

--  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;

if Float_Type'Digits < 20 then return Guess; end if;

--  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;

if Float_Type'Digits < 62 then return Guess; end if;

--  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;

return Guess;

end P_Sum;

end Generic_Pythagorean_Sums;
```