Re: Calculating SQRT in ADA

["David C. Hoos, Sr." <david.c.hoos.sr@ada95.com>]

Marin David Condic wrote in message <36FBCEE2.5DCA764F@pwfl.com>...
>cmcrae wrote:
>>
>>   I am trying to calculate the SQRT in ADA.  Is there a function that I
>> can access?
>> --
>> Posted via Talkway - http://www.talkway.com
>> Surf Usenet at home, on the road, and by email -- always at Talkway.
>
>I certainly hope that the degeneration of this thread into a detailed
>discussion of how to compute SQRT has not discouraged you from using Ada
>:-) I hope you ultimately found a satisfactory answer. (Either the
>numerics appendix in A.5 for Ada95 or compiler specific - post the brand
>and platform & maybe we can help.)
>
>I sometimes wonder if questions like this aren't just tossed out to see
>how many people will run off on a tangent? :-) "Hey everybody! I'm
>looking for the best sorting algorithm available in Ada. Can somebody
>tell me what is the best algorithm?" :-))
>
At the risk of continuing on a tangent ... although the point for Sqrt _per se_ is moot, discussions of quadratic convergence, etc.,
is pertinent to many problems for which the standard libraries don't provide solutions.

To hopefully wrap this up, Ada95 provides a mechanism for cheaply generating the first estimate for square root, which had to be
done in a machine-specific way in Ada83.  The new Ada95 attributes 'Exponent and 'Compose provide the means, as demonstrated in the
example code which follows.

This implementation takes the best of all of the "roll your own" solutions discussed in this thread, and provides a solution in less
than six iterations.

It should be pointed out, however, that the gnat implementation which imports the C library's "sqrt" function is about ten times
faster than this one, on Solaris 7, x86

with Ada.Text_Io;
procedure Test_Sqrt is
   generic
      type Real is digits <>;
   function Sqrt (X : real) return Real;
   function Sqrt (X : real) return real is
      Previous_Previous_Estimate : Real :=
        Real'Compose (0.5, (Real'Exponent (X) + 1) / 2);
      Previous_Estimate          : Real;
      New_Estimate               : Real;
   begin
      if X < 0.0 then
         raise Constraint_Error;
      end if;

      -- Special case Sqrt (0.0) to preserve possible minus sign per
      -- IEEE?? This seeems strange -- why would we ever return a
      -- negative result for a real square root?
      if X = 0.0 then
         return X;
      end if;

      Previous_Estimate := 0.5 *
        (X / Previous_Previous_Estimate + Previous_Previous_Estimate);
      New_Estimate := 0.5 *
        (X / Previous_Estimate + Previous_Estimate);

      while  New_Estimate /= Previous_Estimate and then
        New_Estimate /= Previous_Previous_Estimate loop
         Previous_Previous_Estimate := Previous_Estimate;
         Previous_Estimate := New_Estimate;
         New_Estimate :=
           0.5 * (X / Previous_Estimate + Previous_Estimate);
      end loop;

      return New_Estimate;

   end Sqrt;

   function Long_Float_Sqrt is new Sqrt (Long_Float);

   Root : Long_Float;
   begin
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (2.0)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (0.5)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (2.0e100)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (0.5e-100)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (4.0)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (0.25)));
   end Test_Sqrt;