From: Geert Bosch Date: Thu, 13 Oct 2011 10:51:39 +0000 (+0000) Subject: a-ngrear.adb, [...] (Sqrt): Make generic and move to System.Generic_Array_Operations. X-Git-Tag: releases/gcc-4.7.0~3160 X-Git-Url: http://git.ipfire.org/gitweb.cgi?a=commitdiff_plain;h=edcf5983b69a21cdc76fc31ffbb82598441c6ba3;p=thirdparty%2Fgcc.git a-ngrear.adb, [...] (Sqrt): Make generic and move to System.Generic_Array_Operations. 2011-10-13 Geert Bosch * a-ngrear.adb, s-gearop.adb, s-gearop.ads (Sqrt): Make generic and move to System.Generic_Array_Operations. From-SVN: r179909 --- diff --git a/gcc/ada/ChangeLog b/gcc/ada/ChangeLog index 39d4ec08fbb8..24fd5821f83d 100644 --- a/gcc/ada/ChangeLog +++ b/gcc/ada/ChangeLog @@ -1,3 +1,8 @@ +2011-10-13 Geert Bosch + + * a-ngrear.adb, s-gearop.adb, s-gearop.ads (Sqrt): Make generic and + move to System.Generic_Array_Operations. + 2011-10-13 Geert Bosch * a-ngrear.adb ("abs"): Adjust for modified L2_Norm generic diff --git a/gcc/ada/a-ngrear.adb b/gcc/ada/a-ngrear.adb index 8ffd95e875e3..85c949eebb90 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 -- @@ -119,51 +119,6 @@ package body Ada.Numerics.Generic_Real_Arrays is Y := Old_Y + Sin * (Old_X - Old_Y * Tau); 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 7582e9860e6a..1380cd449cda 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 ca6b7f3586f2..c8eea4f94401 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 -- -----------------