From patchwork Thu Oct 13 10:51:59 2011 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Arnaud Charlet X-Patchwork-Id: 119416 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) by ozlabs.org (Postfix) with SMTP id E6688B6F64 for ; Thu, 13 Oct 2011 21:52:18 +1100 (EST) Received: (qmail 4711 invoked by alias); 13 Oct 2011 10:52:15 -0000 Received: (qmail 4689 invoked by uid 22791); 13 Oct 2011 10:52:14 -0000 X-SWARE-Spam-Status: No, hits=-1.8 required=5.0 tests=AWL,BAYES_00 X-Spam-Check-By: sourceware.org Received: from rock.gnat.com (HELO rock.gnat.com) (205.232.38.15) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Thu, 13 Oct 2011 10:52:00 +0000 Received: from localhost (localhost.localdomain [127.0.0.1]) by filtered-rock.gnat.com (Postfix) with ESMTP id 3CBA52BB4BA; Thu, 13 Oct 2011 06:51:59 -0400 (EDT) Received: from rock.gnat.com ([127.0.0.1]) by localhost (rock.gnat.com [127.0.0.1]) (amavisd-new, port 10024) with LMTP id 1G3dXq0Sqb8W; Thu, 13 Oct 2011 06:51:59 -0400 (EDT) Received: from kwai.gnat.com (kwai.gnat.com [205.232.38.4]) by rock.gnat.com (Postfix) with ESMTP id 29AA92BB4B8; Thu, 13 Oct 2011 06:51:59 -0400 (EDT) Received: by kwai.gnat.com (Postfix, from userid 4192) id 2897F92BF6; Thu, 13 Oct 2011 06:51:59 -0400 (EDT) Date: Thu, 13 Oct 2011 06:51:59 -0400 From: Arnaud Charlet To: gcc-patches@gcc.gnu.org Cc: Geert Bosch Subject: [Ada] Make local Sqrt implementation generic Message-ID: <20111013105159.GA4837@adacore.com> MIME-Version: 1.0 Content-Disposition: inline User-Agent: Mutt/1.5.20 (2009-06-14) Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Delivered-To: mailing list gcc-patches@gcc.gnu.org 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 * a-ngrear.adb, s-gearop.adb, s-gearop.ads (Sqrt): Make generic and move to System.Generic_Array_Operations. 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 -- -----------------