Integer Square Root

Integer Square Root

[Matthew Heaney]

W. Wesley Groleau <mailto:wwgrol@pseserv3.fw.hac.com> pointed out to me
some errors in an algorithm I posted for performing the square root of an
integer number without using floating point arithmetic.

That I did that off the top of my head, without any testing (neither static
nor dynamic), and the fact that that algorithm basically doesn't work,
should prove to everyone (especially me) that "off-the-cuff" programming
doesn't work!  Thank you, Wes, for pointing out to me the errors of my
ways.

So here, if I may be given a chance to right my wrongs, is a "correct" version:

   function Square_Root (N : Natural) return Natural is
      X0 : Natural := N / 2;
      X1 : Natural := 2;
   begin
      case N is 
         when 0 =>
            return 0;

         when 1 =>
           return 1;

         when others =>
             while abs (X0 - X1) > 1 loop
                  X0 := (X0 + X1) / 2;
                  X1 := N / X0;
            end loop;

            return X0;

      end case;
   end Square_Root;

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.

matt

--------------------------------------------------------------------
Matthew Heaney
Software Development Consultant
mheaney@ni.net
(818) 985-1271

Re: Integer Square Root

[Keith Thompson]

In <mheaney-ya023180001510962133390001@news.ni.net> mheaney@ni.net (Matthew Heaney) writes:
[...]
> That I did that off the top of my head, without any testing (neither static
> nor dynamic), and the fact that that algorithm basically doesn't work,
> should prove to everyone (especially me) that "off-the-cuff" programming
> doesn't work!  Thank you, Wes, for pointing out to me the errors of my
> ways.
> 
> So here, if I may be given a chance to right my wrongs, is a "correct" version:
> 
>    function Square_Root (N : Natural) return Natural is
[...]

Sorry, but this still doesn't work.  It says the square root of 6 is 3.
If it's supposed to truncate, Square_Root(6) should return 2; if it's
supposed to round, sqrt(6.0) = 2.4495 (approximately), so Square_Root(6)
should still return 2.  It also doesn't always return increasing values
for increasing numbers; Square_Root(135) returns 12 and Square_Root(136)
returns 11.

-- 
Keith Thompson (The_Other_Keith) kst@thomsoft.com <*>
TeleSoft^H^H^H^H^H^H^H^H Alsys^H^H^H^H^H Thomson Software Products
10251 Vista Sorrento Parkway, Suite 300, San Diego, CA, USA, 92121-2706
FIJAGDWOL

Re: Integer Square Root

[Matthew Heaney]

In article <DzEEq3.9tI@thomsoft.com>, kst@thomsoft.com (Keith Thompson) wrote:

>>    function Square_Root (N : Natural) return Natural is
>[...]
>
>Sorry, but this still doesn't work.  It says the square root of 6 is 3.
>If it's supposed to truncate, Square_Root(6) should return 2; if it's
>supposed to round, sqrt(6.0) = 2.4495 (approximately), so Square_Root(6)
>should still return 2.  It also doesn't always return increasing values
>for increasing numbers; Square_Root(135) returns 12 and Square_Root(136)
>returns 11.

Well, yes and no.  I intentionally wrote the return statement in the others
branch to be a "don't care"; it wasn't "supposed" to be truncation or
rounding.   However, there are obvious virtues to having a continuously
increasing square root function!

If you want truncation behavior, the modification is simple:

   return Natural'Min (X0, X1);

For rounding behavior, 

   return Natural'Max (X0, X1);

Perhaps I should have stated the behavior explicitly; thank you for
pointing this out.

--------------------------------------------------------------------
Matthew Heaney
Software Development Consultant
mheaney@ni.net
(818) 985-1271

Re: Integer Square Root

[Oliver Kellogg]

-- This is a square root algorithm for unsigned 32 bit integers that
-- uses addition, subtraction, and shifts only.
-- Result is *truncated*, i.e. Sqrt(62) = 7.
-- For a description of the algorithm, see the German computer magazine
-- c't, January 1990, page 300 ff. (Otto Peter, "Prozessor zieht Wurzeln")

with Interfaces;
use  Interfaces;
function Sqrt (X : Unsigned_32) return Unsigned_32 is
   Root : Unsigned_32 := 0;
   M    : Unsigned_32 := 16#4000_0000#;
   X1   : Unsigned_32 := X;
   X2   : Unsigned_32;
begin
   loop
      X2 := Root + M;
      Root := Shift_Right (Root, 1);
      if X2 <= X1 then
         X1 := X1 - X2;
         Root := Root + M;
      end if;
      M := Shift_Right (M, 2);
      exit when M = 0;
   end loop;
   return Root;
end Sqrt;

Re: Integer Square Root

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

Integer Square Root

[W. Wesley Groleau (Wes)]

I meant to post this a long time ago.  Sorry about the delay.

M. Heaney's integer square root algorithm should not be used without testing.
After applying the correction he posted, here's how it performs:

If the input is zero or one, a divide by zero exception is raised.
(constraint_error)

perfect square:
If the input is N*N where N is any integer greater than one, the answer
is correct (after applying the correction you posted).

perfect square minus one:
For the same N, an input of N*N-1 never terminates.

For all other positive integers, the answer is equivalent to the float
answer with the fractional part truncated, i.e., root 62 comes out 7 where
integer(sqrt(float(62))) would be 8.  (Some people may prefer truncation).

For negative numbers, sometimes it produces an answer (possibly a negative
answer), sometimes it doesn't terminate, sometimes it raises the exception.

---------------------------------------------------------------------------
W. Wesley Groleau (Wes)                                Office: 219-429-4923
Hughes Defense Communications (MS 10-40)                 Home: 219-471-7206
Fort Wayne,  IN   46808                  (Unix): wwgrol@pseserv3.fw.hac.com
---------------------------------------------------------------------------