--  (C) Copyright 1999 by John Halleck
--  All rights reserved

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

   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;