Patchwork [Ada] Make local Sqrt implementation generic

login
register
mail settings
Submitter Arnaud Charlet
Date Oct. 13, 2011, 10:51 a.m.
Message ID <20111013105159.GA4837@adacore.com>
Download mbox | patch
Permalink /patch/119416/
State New
Headers show

Comments

Arnaud Charlet - Oct. 13, 2011, 10:51 a.m.
This prepares for reusing the Sqrt implementation in Generic_Complex_Arrays.
The local implementation avoids having to instantiate entire new copies of
Generic_Elementary_Functions just to get square root.

Tested on x86_64-pc-linux-gnu, committed on trunk

2011-10-13  Geert Bosch  <bosch@adacore.com>

	* a-ngrear.adb, s-gearop.adb, s-gearop.ads (Sqrt): Make generic and
	move to System.Generic_Array_Operations.

Patch

Index: a-ngrear.adb
===================================================================
--- a-ngrear.adb	(revision 179908)
+++ a-ngrear.adb	(working copy)
@@ -102,10 +102,10 @@ 
    procedure Swap (Left, Right : in out Real);
    --  Exchange Left and Right
 
-   function Sqrt (X : Real) return Real;
-   --  Sqrt is implemented locally here, in order to avoid dragging in all of
-   --  the elementary functions. Speed of the square root is not a big concern
-   --  here. This also avoids depending on a specific floating point type.
+   function Sqrt is new Ops.Sqrt (Real);
+   --  Instant a generic square root implementation here, in order to avoid
+   --  instantiating a complete copy of Generic_Elementary_Functions.
+   --  Speed of the square root is not a big concern here.
 
    ------------
    -- Rotate --
@@ -120,51 +120,6 @@ 
    end Rotate;
 
    ----------
-   -- Sqrt --
-   ----------
-
-   function Sqrt (X : Real) return Real is
-      Root, Next : Real;
-
-   begin
-      --  Be defensive: any comparisons with NaN values will yield False.
-
-      if not (X > 0.0) then
-         if X = 0.0 then
-            return X;
-         else
-            raise Argument_Error;
-         end if;
-      end if;
-
-      --  Compute an initial estimate based on:
-
-      --     X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
-
-      --  where M is the mantissa, R is the radix and E the exponent.
-
-      --  By ignoring the mantissa and ignoring the case of an odd
-      --  exponent, we get a final error that is at most R. In other words,
-      --  the result has about a single bit precision.
-
-      Root := Real (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
-
-      --  Because of the poor initial estimate, use the Babylonian method of
-      --  computing the square root, as it is stable for all inputs. Every step
-      --  will roughly double the precision of the result. Just a few steps
-      --  suffice in most cases. Eight iterations should give about 2**8 bits
-      --  of precision.
-
-      for J in 1 .. 8 loop
-         Next := (Root + X / Root) / 2.0;
-         exit when Root = Next;
-         Root := Next;
-      end loop;
-
-      return Root;
-   end Sqrt;
-
-   ----------
    -- Swap --
    ----------
 
Index: s-gearop.adb
===================================================================
--- s-gearop.adb	(revision 179908)
+++ s-gearop.adb	(working copy)
@@ -29,6 +29,8 @@ 
 --                                                                          --
 ------------------------------------------------------------------------------
 
+with Ada.Numerics; use Ada.Numerics;
+
 package body System.Generic_Array_Operations is
 
    --  The local function Check_Unit_Last computes the index
@@ -567,6 +569,56 @@ 
       return R;
    end Scalar_Vector_Elementwise_Operation;
 
+   ----------
+   -- Sqrt --
+   ----------
+
+   function Sqrt (X : Real'Base) return Real'Base is
+      Root, Next : Real'Base;
+
+   begin
+      --  Be defensive: any comparisons with NaN values will yield False.
+
+      if not (X > 0.0) then
+         if X = 0.0 then
+            return X;
+         else
+            raise Argument_Error;
+         end if;
+
+      elsif X > Real'Base'Last then
+         --  X is infinity, which is its own square root
+
+         return X;
+      end if;
+
+      --  Compute an initial estimate based on:
+
+      --     X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
+
+      --  where M is the mantissa, R is the radix and E the exponent.
+
+      --  By ignoring the mantissa and ignoring the case of an odd
+      --  exponent, we get a final error that is at most R. In other words,
+      --  the result has about a single bit precision.
+
+      Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
+
+      --  Because of the poor initial estimate, use the Babylonian method of
+      --  computing the square root, as it is stable for all inputs. Every step
+      --  will roughly double the precision of the result. Just a few steps
+      --  suffice in most cases. Eight iterations should give about 2**8 bits
+      --  of precision.
+
+      for J in 1 .. 8 loop
+         Next := (Root + X / Root) / 2.0;
+         exit when Root = Next;
+         Root := Next;
+      end loop;
+
+      return Root;
+   end Sqrt;
+
    ---------------------------
    -- Matrix_Matrix_Product --
    ---------------------------
Index: s-gearop.ads
===================================================================
--- s-gearop.ads	(revision 179908)
+++ s-gearop.ads	(working copy)
@@ -388,6 +388,14 @@ 
      (Left  : Left_Matrix;
       Right : Right_Matrix) return Result_Matrix;
 
+   ----------
+   -- Sqrt --
+   ----------
+
+   generic
+      type Real is digits <>;
+   function Sqrt (X : Real'Base) return Real'Base;
+
    -----------------
    -- Swap_Column --
    -----------------