From patchwork Tue Jun 22 17:12:07 2010 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Arnaud Charlet X-Patchwork-Id: 56538 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 18D75B6F06 for ; Wed, 23 Jun 2010 03:12:16 +1000 (EST) Received: (qmail 6499 invoked by alias); 22 Jun 2010 17:12:15 -0000 Received: (qmail 6489 invoked by uid 22791); 22 Jun 2010 17:12:14 -0000 X-SWARE-Spam-Status: No, hits=-1.8 required=5.0 tests=AWL, BAYES_00, T_RP_MATCHES_RCVD X-Spam-Check-By: sourceware.org Received: from mel.act-europe.fr (HELO mel.act-europe.fr) (212.99.106.210) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 22 Jun 2010 17:12:05 +0000 Received: from localhost (localhost [127.0.0.1]) by filtered-smtp.eu.adacore.com (Postfix) with ESMTP id 4158FCB026C; Tue, 22 Jun 2010 19:12:07 +0200 (CEST) Received: from mel.act-europe.fr ([127.0.0.1]) by localhost (smtp.eu.adacore.com [127.0.0.1]) (amavisd-new, port 10024) with ESMTP id 47sQQpJ6nZob; Tue, 22 Jun 2010 19:12:07 +0200 (CEST) Received: from saumur.act-europe.fr (saumur.act-europe.fr [10.10.0.183]) by mel.act-europe.fr (Postfix) with ESMTP id 2EC75CB023D; Tue, 22 Jun 2010 19:12:07 +0200 (CEST) Received: by saumur.act-europe.fr (Postfix, from userid 525) id 2727BD9BB4; Tue, 22 Jun 2010 19:12:07 +0200 (CEST) Date: Tue, 22 Jun 2010 19:12:07 +0200 From: Arnaud Charlet To: gcc-patches@gcc.gnu.org Cc: Geert Bosch Subject: [Ada] Fix scaled complex multiplication in case of overflow Message-ID: <20100622171207.GA15102@adacore.com> Mime-Version: 1.0 Content-Disposition: inline User-Agent: Mutt/1.5.9i X-IsSubscribed: yes 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 In case the canonical implementation of complex multiplication would overflow, an attempt was made to use scaling. However, the scaling would not be done in case of NaN results that may result from cancelling infinities and the scaling would not be sufficient in most cases. Also scaling for the real part used a wrong sign. This patch addresses all issues. -- The following test case must compile and execute quietly. with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types; procedure Huge is L : constant Float := Float'Pred (Float'Pred (Float'Pred (Float'Last))); -- Largest number with 2 trailing bits equal to 0 X1 : constant Complex := (L, L); Y1 : constant Complex := (3.0, 2.0); X2 : constant Complex := (L, -2.0); Y2 : constant Complex := (L, 3.0); X1Y1 : constant Float := Re (X1 * Y1); X2Y2 : constant Float := Im (X2 * Y2); begin if X1Y1 /= L or else X2Y2 /= L then raise Program_Error; end if; end Huge; Tested on x86_64-pc-linux-gnu, committed on trunk 2010-06-22 Geert Bosch * a-ngcoty.adb ("*"): Rewrite complex multiplication to use proper scaling in case of overflow or NaN results. Index: a-ngcoty.adb =================================================================== --- a-ngcoty.adb (revision 161073) +++ a-ngcoty.adb (working copy) @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 1992-2009, Free Software Foundation, Inc. -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- @@ -43,6 +43,12 @@ package body Ada.Numerics.Generic_Comple --------- function "*" (Left, Right : Complex) return Complex is + + Scale : constant R := R (R'Machine_Radix) ** ((R'Machine_Emax - 1) / 2); + -- In case of overflow, scale the operands by the largest power of the + -- radix (to avoid rounding error), so that the square of the scale does + -- not overflow itself. + X : R; Y : R; @@ -53,14 +59,19 @@ package body Ada.Numerics.Generic_Comple -- If either component overflows, try to scale (skip in fast math mode) if not Standard'Fast_Math then - if abs (X) > R'Last then - X := R'(4.0) * (R'(Left.Re / 2.0) * R'(Right.Re / 2.0) - - R'(Left.Im / 2.0) * R'(Right.Im / 2.0)); + + -- ??? the test below is weird, it needs a comment, otherwise I or + -- someone else will change it back to R'Last > abs (X) ??? + + if not (abs (X) <= R'Last) then + X := Scale**2 * ((Left.Re / Scale) * (Right.Re / Scale) - + (Left.Im / Scale) * (Right.Im / Scale)); end if; - if abs (Y) > R'Last then - Y := R'(4.0) * (R'(Left.Re / 2.0) * R'(Right.Im / 2.0) - - R'(Left.Im / 2.0) * R'(Right.Re / 2.0)); + -- ??? same weird test ??? + if not (abs (Y) <= R'Last) then + Y := Scale**2 * ((Left.Re / Scale) * (Right.Im / Scale) + + (Left.Im / Scale) * (Right.Re / Scale)); end if; end if; @@ -569,7 +580,8 @@ package body Ada.Numerics.Generic_Comple -- in order to prevent inaccuracies on machines where not all -- immediate expressions are rounded, such as PowerPC. - if Re2 > R'Last then + -- ??? same weird test, why not Re2 > R'Last ??? + if not (Re2 <= R'Last) then raise Constraint_Error; end if; @@ -582,7 +594,8 @@ package body Ada.Numerics.Generic_Comple begin Im2 := X.Im ** 2; - if Im2 > R'Last then + -- ??? same weird test + if not (Im2 <= R'Last) then raise Constraint_Error; end if;