-- --
-- B o d y --
-- --
--- Copyright (C) 1992-2005 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- --
--- ware Foundation; either version 2, or (at your option) any later ver- --
+-- ware Foundation; either version 3, or (at your option) any later ver- --
-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
--- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
--- for more details. You should have received a copy of the GNU General --
--- Public License distributed with GNAT; see file COPYING. If not, write --
--- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
--- MA 02111-1307, USA. --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
-- --
--- As a special exception, if other files instantiate generics from this --
--- unit, or you link this unit with other files to produce an executable, --
--- this unit does not by itself cause the resulting executable to be --
--- covered by the GNU General Public License. This exception does not --
--- however invalidate any other reasons why the executable file might be --
--- covered by the GNU Public License. --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
+-- --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
-- --
-- GNAT was originally developed by the GNAT team at New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc. --
------------------------------------------------------------------------------
with Ada.Numerics.Aux; use Ada.Numerics.Aux;
+
package body Ada.Numerics.Generic_Complex_Types is
subtype R is Real'Base;
---------
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;
X := Left.Re * Right.Re - Left.Im * Right.Im;
Y := Left.Re * Right.Im + Left.Im * Right.Re;
- -- If either component overflows, try to scale.
+ -- If either component overflows, try to scale (skip in fast math mode)
- 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));
- end if;
+ if not Standard'Fast_Math then
+
+ -- Note that the test below is written as a negation. This is to
+ -- account for the fact that X and Y may be NaNs, because both of
+ -- their operands could overflow. Given that all operations on NaNs
+ -- return false, the test can only be written thus.
- 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));
+ 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 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;
return (X, Y);
function "*" (Left, Right : Imaginary) return Real'Base is
begin
- return -R (Left) * R (Right);
+ return -(R (Left) * R (Right));
end "*";
function "*" (Left : Complex; Right : Real'Base) return Complex is
-- 1.0 / infinity, and the closest model number will be zero.
begin
-
while Exp /= 0 loop
if Exp rem 2 /= 0 then
Result := Result * Factor;
return R'(1.0) / Result;
exception
-
when Constraint_Error =>
return (0.0, 0.0);
end;
c : constant R := Right.Re;
d : constant R := Right.Im;
begin
- return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
- Im => -(a * d) / (c ** 2 + d ** 2));
+ return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
+ Im => -((a * d) / (c ** 2 + d ** 2)));
end "/";
function "/" (Left : Complex; Right : Imaginary) return Complex is
d : constant R := R (Right);
begin
- return (b / d, -a / d);
+ return (b / d, -(a / d));
end "/";
function "/" (Left : Imaginary; Right : Complex) return Complex is
function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
begin
- return Imaginary (-Left / R (Right));
+ return Imaginary (-(Left / R (Right)));
end "/";
---------
-- 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;
begin
Im2 := X.Im ** 2;
- if Im2 > R'Last then
+ -- ??? same weird test
+ if not (Im2 <= R'Last) then
raise Constraint_Error;
end if;
elsif Im2 = 0.0 then
return abs (X.Re);
- -- in all other cases, the naive computation will do.
+ -- In all other cases, the naive computation will do
else
return R (Sqrt (Double (Re2 + Im2)));
-- Set_Im --
------------
- procedure Set_Im (X : in out Complex; Im : in Real'Base) is
+ procedure Set_Im (X : in out Complex; Im : Real'Base) is
begin
X.Im := Im;
end Set_Im;
- procedure Set_Im (X : out Imaginary; Im : in Real'Base) is
+ procedure Set_Im (X : out Imaginary; Im : Real'Base) is
begin
X := Imaginary (Im);
end Set_Im;
-- Set_Re --
------------
- procedure Set_Re (X : in out Complex; Re : in Real'Base) is
+ procedure Set_Re (X : in out Complex; Re : Real'Base) is
begin
X.Re := Re;
end Set_Re;