------------------------------------------------------------------------------
-- --
--- GNAT RUNTIME COMPONENTS --
+-- GNAT RUN-TIME COMPONENTS --
-- --
-- A D A . N U M E R I C S . G E N E R I C _ C O M P L E X _ T Y P E S --
-- --
-- B o d y --
-- --
--- $Revision: 1.16 $
--- --
--- Copyright (C) 1992-2001 Free Software Foundation, Inc. --
+-- Copyright (C) 1992-2006, 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- --
-- 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. --
+-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
+-- Boston, MA 02110-1301, USA. --
-- --
-- As a special exception, if other files instantiate generics from this --
-- unit, or you link this unit with other files to produce an executable, --
-- covered by the GNU Public License. --
-- --
-- GNAT was originally developed by the GNAT team at New York University. --
--- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
+-- 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;
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
if abs (X) > R'Last then
- X := R' (4.0) * (R'(Left.Re / 2.0) * R'(Right.Re / 2.0)
+ 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 abs (Y) > R'Last then
- Y := R' (4.0) * (R'(Left.Re / 2.0) * R'(Right.Im / 2.0)
+ Y := R'(4.0) * (R'(Left.Re / 2.0) * R'(Right.Im / 2.0)
- R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
end if;
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
Exp := Exp / 2;
end loop;
- return R ' (1.0) / Result;
+ return R'(1.0) / Result;
exception
end "**";
function "**" (Left : Imaginary; Right : Integer) return Complex is
- M : R := R (Left) ** Right;
+ M : constant R := R (Left) ** Right;
begin
case Right mod 4 is
when 0 => return (M, 0.0);
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 "/";
---------
-- we can use an explicit comparison to determine whether to use
-- the scaling expression.
+ -- The scaling expression is computed in double format throughout
+ -- in order to prevent inaccuracies on machines where not all
+ -- immediate expressions are rounded, such as PowerPC.
+
if Re2 > R'Last then
raise Constraint_Error;
end if;
exception
when Constraint_Error =>
- return abs (X.Re)
- * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
+ return R (Double (abs (X.Re))
+ * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
end;
begin
exception
when Constraint_Error =>
- return abs (X.Im)
- * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
+ return R (Double (abs (X.Im))
+ * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
end;
-- Now deal with cases of underflow. If only one of the squares
else
if abs (X.Re) > abs (X.Im) then
return
- abs (X.Re)
- * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
+ R (Double (abs (X.Re))
+ * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
else
return
- abs (X.Im)
- * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
+ R (Double (abs (X.Im))
+ * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
end if;
end if;
return abs (X.Im);
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;
OSDN Git Service
