Re: Integer Square Root

[Mats.Weber@elca-matrix.ch (Mats Weber)]

>This is similar to an algorithm for fixed point types that appears in
>Section 5.4 of the Ada 83 Rationale.
>
>If anyone is interested in a generic version that works for any integer
>type, then let me know, and I'll post it.

Here is the one I did a long time ago. I have tested it quite extensively
and it works (always finishes in finite time).

-- SQUARE ROOT FUNCTIONS
   ---------------------

-- Creation : 17-NOV-1989 by Mats Weber, taken from package Math_Pack.


package Square_Root_Functions is
-----------------------------

   generic
      type Int is range <>;
   function Integer_Sqrt (X : Int) return Int;

   generic
      type Real is delta <>;
   function Fixed_Sqrt (X : Real) return Real;

   Sqrt_Of_Negative_Value : exception;

end Square_Root_Functions;


-- SQUARE ROOT FUNCTIONS
   ---------------------

-- Creation : 17-NOV-1989 by Mats Weber.


package body Square_Root_Functions is
----------------------------------

   function Integer_Sqrt (X : Int) return Int is

      U     : Int := X;
      New_U : Int;

   begin
      if X < 0 then
         raise Sqrt_Of_Negative_Value;
      end if;
      if X = 0 then
         return 0;
      end if;
      loop
         New_U := (U + X/U)/2;
         exit when New_U >= U;
         U := New_U;
      end loop;
      return U;
   end Integer_Sqrt;


   function Fixed_Sqrt (X : Real) return Real is

      U     : Real := X;
      New_U : Real;

   begin
      if X < 0.0 then
         raise Sqrt_Of_Negative_Value;
      end if;
      if X = 0.0 then
         return 0.0;
      end if;
      if X >= 1.0 then          -- the series is decreasing
         loop
            New_U := (U + Real(X/U))/2;
            exit when New_U >= U;
            U := New_U;
         end loop;
         return U;
      else                      -- the series is increasing
         loop
            New_U := (U + Real(X/U))/2;
            exit when New_U <= U;
            U := New_U;
         end loop;
         return New_U;
      end if;
   end Fixed_Sqrt;

end Square_Root_Functions;