Re: Calculating SQRT in ADA

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

Hans Marqvardsen wrote in message <36F8D9B7.1F8D@ddre.dk>...
>John Herro wrote:
>>
>> Robert Dewar wrote:
>> >>   while ...
>> >>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
>> > Why the epsilon test, this iteration will terminate to
>> > exact equality in any reasonable floating-point model ...
>>
>> I don't know if the Open Ada 83 compiler has a "reasonable" floating-point
>> model, but I do know that if I test for exact equality, instantiate with type
>> Float, and compute Sqrt(2.0) with that compiler, the program goes into an
>> infinite loop.  I think that in general, you can't guarantee that Guess*Guess
>> will ever exactly equal X.
>>
>
>Clearly, in any floating-point model the square root of X can only
>be represented  exactly  for very few values of X.
>
>However any "reasonable" model may be expected to contain two
>neighbour numbers G and H such that
>
>   G*G <= X <= H*H
>
>Then one may expect convergence to either
>G, or H, or to the alternating sequence G,H,G,H,...
>
>Hence you need the epsilon test after all,
>but probably the constant 3.0 can safely be reduced.

I take issue with all of the previous answers on at least one point in each.

First, the test should not be for Guess * Guess / X = 1.0, because this is far too expensive of CPU cycles and unnecessary.  The
correct test, as Hans points out is for equality of successive estimates, or alternating.  This permits the solution to be obtained
with one division, one multiplication, one addition, and one or two floating point comparisons per iteration.

So, no epsilon test _per se_ is required.

As Robert correctly points out the choice of the initial estimate is important.  One solution to this is to examine the exponent of
the value of which the root is desired (argument), and start with an estimate of the root having an exponent half that of the
argument.

For example, of the machine radix is 2, then if the argument has an odd exponent, one can denormalize the number by increasing the
exponent by one (so it can now be evenly divisible by 2), and halving the mantissa.
This puts the mantissa of the denormalized argument on the range 0.25 .. 1.0, so its square root will be the correct value for the
mantissa of the square root, and the exponent of the square root is simply half of the (possibly incremented) exponent of the
argument.

This approach, however, is machine-dependent, but we actually found this approach faster than any other, on a machine whose
OS-provided math library (which was used by all language compilers) provided abysmal performance.

The square root of the mantissa can be closely estimated using a Pade' approximation, such that generally no more than two
Newton-Raphson iterations are required for convergence.

I became very concerned about this sort of thing about thirty years ago when CPU time cost me 40 cents per second on a GE 635.

Here is an implementation which is machine-independent, but suffers from the initial estimate problem.  It will, however,
unconditionally converge.


   generic
      type Real is digits <>;
   function Sqrt (X : real) return Real;

   function Sqrt (X : real) return real is
      Previous_Estimate : Real := X;
      New_Estimate : Real;
      Increasing : Boolean;
      Reversed : Boolean := False;
      Alternating : Boolean := False;
   begin
      if X < 0.0 then
         raise Constraint_Error;
      end if;
      if X = 0.0 then
         return 0.0;
      end if;
      New_Estimate :=
        0.5 * (X / Previous_Estimate + Previous_Estimate);
      Increasing := New_Estimate > Previous_Estimate;
      while  not alternating and then
        New_Estimate /= Previous_estimate loop
         if not Reversed then
            Reversed :=
              Increasing /= (New_Estimate > Previous_Estimate);
         else
            Alternating :=
              Increasing = (New_Estimate > Previous_Estimate);
         end if;
         Previous_Estimate := New_Estimate;
         New_Estimate :=
           0.5 * (X / Previous_Estimate + Previous_Estimate);
      end loop;
      return New_Estimate;
   end Sqrt;