a-ngrear.adb, [...] (Sqrt): Make generic and move to System.Generic_Array_Operations.

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.

From-SVN: r179909

4 files changed, 69 insertions, 49 deletions

diff --git a/gcc/ada/ChangeLog b/gcc/ada/ChangeLog
index 39d4ec0..24fd582 100644
--- a/gcc/ada/ChangeLog
+++ b/gcc/ada/ChangeLog
@@ -1,5 +1,10 @@
 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.
+
+2011-10-13  Geert Bosch  <bosch@adacore.com>
+
 	* a-ngrear.adb ("abs"): Adjust for modified L2_Norm generic
 	* s-gearop.ads (L2_Norm): Change profile to be suitable for
 	Complex_Vector
diff --git a/gcc/ada/a-ngrear.adb b/gcc/ada/a-ngrear.adb
index 8ffd95e..85c949e 100644
--- a/gcc/ada/a-ngrear.adb
+++ b/gcc/ada/a-ngrear.adb
@@ -102,10 +102,10 @@ package body Ada.Numerics.Generic_Real_Arrays is
    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 @@ package body Ada.Numerics.Generic_Real_Arrays is
    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 --
    ----------
 
diff --git a/gcc/ada/s-gearop.adb b/gcc/ada/s-gearop.adb
index 7582e98..1380cd4 100644
--- a/gcc/ada/s-gearop.adb
+++ b/gcc/ada/s-gearop.adb
@@ -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 @@ package body System.Generic_Array_Operations is
       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 --
    ---------------------------
diff --git a/gcc/ada/s-gearop.ads b/gcc/ada/s-gearop.ads
index ca6b7f3..c8eea4f 100644
--- a/gcc/ada/s-gearop.ads
+++ b/gcc/ada/s-gearop.ads
@@ -388,6 +388,14 @@ pragma Pure (Generic_Array_Operations);
      (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 --
    -----------------