diff mbox

[Ada] Fix scaled complex multiplication in case of overflow

Message ID 20100622171207.GA15102@adacore.com
State New
Headers show

Commit Message

Arnaud Charlet June 22, 2010, 5:12 p.m. UTC 
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  <bosch@adacore.com>

	* a-ngcoty.adb ("*"): Rewrite complex multiplication to use proper
	scaling in case of overflow or NaN results.
diff mbox

Patch

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;