682 lines
19 KiB
Ada
682 lines
19 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- 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 --
|
|
-- --
|
|
-- 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 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. --
|
|
-- --
|
|
-- 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;
|
|
|
|
Two_Pi : constant R := R (2.0) * Pi;
|
|
Half_Pi : constant R := Pi / R (2.0);
|
|
|
|
---------
|
|
-- "*" --
|
|
---------
|
|
|
|
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;
|
|
|
|
begin
|
|
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 (skip in fast math mode)
|
|
|
|
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 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);
|
|
end "*";
|
|
|
|
function "*" (Left, Right : Imaginary) return Real'Base is
|
|
begin
|
|
return -(R (Left) * R (Right));
|
|
end "*";
|
|
|
|
function "*" (Left : Complex; Right : Real'Base) return Complex is
|
|
begin
|
|
return Complex'(Left.Re * Right, Left.Im * Right);
|
|
end "*";
|
|
|
|
function "*" (Left : Real'Base; Right : Complex) return Complex is
|
|
begin
|
|
return (Left * Right.Re, Left * Right.Im);
|
|
end "*";
|
|
|
|
function "*" (Left : Complex; Right : Imaginary) return Complex is
|
|
begin
|
|
return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
|
|
end "*";
|
|
|
|
function "*" (Left : Imaginary; Right : Complex) return Complex is
|
|
begin
|
|
return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
|
|
end "*";
|
|
|
|
function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
|
|
begin
|
|
return Left * Imaginary (Right);
|
|
end "*";
|
|
|
|
function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
|
|
begin
|
|
return Imaginary (Left * R (Right));
|
|
end "*";
|
|
|
|
----------
|
|
-- "**" --
|
|
----------
|
|
|
|
function "**" (Left : Complex; Right : Integer) return Complex is
|
|
Result : Complex := (1.0, 0.0);
|
|
Factor : Complex := Left;
|
|
Exp : Integer := Right;
|
|
|
|
begin
|
|
-- We use the standard logarithmic approach, Exp gets shifted right
|
|
-- testing successive low order bits and Factor is the value of the
|
|
-- base raised to the next power of 2. For positive exponents we
|
|
-- multiply the result by this factor, for negative exponents, we
|
|
-- divide by this factor.
|
|
|
|
if Exp >= 0 then
|
|
|
|
-- For a positive exponent, if we get a constraint error during
|
|
-- this loop, it is an overflow, and the constraint error will
|
|
-- simply be passed on to the caller.
|
|
|
|
while Exp /= 0 loop
|
|
if Exp rem 2 /= 0 then
|
|
Result := Result * Factor;
|
|
end if;
|
|
|
|
Factor := Factor * Factor;
|
|
Exp := Exp / 2;
|
|
end loop;
|
|
|
|
return Result;
|
|
|
|
else -- Exp < 0 then
|
|
|
|
-- For the negative exponent case, a constraint error during this
|
|
-- calculation happens if Factor gets too large, and the proper
|
|
-- response is to return 0.0, since what we essentially have 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;
|
|
end if;
|
|
|
|
Factor := Factor * Factor;
|
|
Exp := Exp / 2;
|
|
end loop;
|
|
|
|
return R'(1.0) / Result;
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
return (0.0, 0.0);
|
|
end;
|
|
end if;
|
|
end "**";
|
|
|
|
function "**" (Left : Imaginary; Right : Integer) return Complex is
|
|
M : constant R := R (Left) ** Right;
|
|
begin
|
|
case Right mod 4 is
|
|
when 0 => return (M, 0.0);
|
|
when 1 => return (0.0, M);
|
|
when 2 => return (-M, 0.0);
|
|
when 3 => return (0.0, -M);
|
|
when others => raise Program_Error;
|
|
end case;
|
|
end "**";
|
|
|
|
---------
|
|
-- "+" --
|
|
---------
|
|
|
|
function "+" (Right : Complex) return Complex is
|
|
begin
|
|
return Right;
|
|
end "+";
|
|
|
|
function "+" (Left, Right : Complex) return Complex is
|
|
begin
|
|
return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
|
|
end "+";
|
|
|
|
function "+" (Right : Imaginary) return Imaginary is
|
|
begin
|
|
return Right;
|
|
end "+";
|
|
|
|
function "+" (Left, Right : Imaginary) return Imaginary is
|
|
begin
|
|
return Imaginary (R (Left) + R (Right));
|
|
end "+";
|
|
|
|
function "+" (Left : Complex; Right : Real'Base) return Complex is
|
|
begin
|
|
return Complex'(Left.Re + Right, Left.Im);
|
|
end "+";
|
|
|
|
function "+" (Left : Real'Base; Right : Complex) return Complex is
|
|
begin
|
|
return Complex'(Left + Right.Re, Right.Im);
|
|
end "+";
|
|
|
|
function "+" (Left : Complex; Right : Imaginary) return Complex is
|
|
begin
|
|
return Complex'(Left.Re, Left.Im + R (Right));
|
|
end "+";
|
|
|
|
function "+" (Left : Imaginary; Right : Complex) return Complex is
|
|
begin
|
|
return Complex'(Right.Re, R (Left) + Right.Im);
|
|
end "+";
|
|
|
|
function "+" (Left : Imaginary; Right : Real'Base) return Complex is
|
|
begin
|
|
return Complex'(Right, R (Left));
|
|
end "+";
|
|
|
|
function "+" (Left : Real'Base; Right : Imaginary) return Complex is
|
|
begin
|
|
return Complex'(Left, R (Right));
|
|
end "+";
|
|
|
|
---------
|
|
-- "-" --
|
|
---------
|
|
|
|
function "-" (Right : Complex) return Complex is
|
|
begin
|
|
return (-Right.Re, -Right.Im);
|
|
end "-";
|
|
|
|
function "-" (Left, Right : Complex) return Complex is
|
|
begin
|
|
return (Left.Re - Right.Re, Left.Im - Right.Im);
|
|
end "-";
|
|
|
|
function "-" (Right : Imaginary) return Imaginary is
|
|
begin
|
|
return Imaginary (-R (Right));
|
|
end "-";
|
|
|
|
function "-" (Left, Right : Imaginary) return Imaginary is
|
|
begin
|
|
return Imaginary (R (Left) - R (Right));
|
|
end "-";
|
|
|
|
function "-" (Left : Complex; Right : Real'Base) return Complex is
|
|
begin
|
|
return Complex'(Left.Re - Right, Left.Im);
|
|
end "-";
|
|
|
|
function "-" (Left : Real'Base; Right : Complex) return Complex is
|
|
begin
|
|
return Complex'(Left - Right.Re, -Right.Im);
|
|
end "-";
|
|
|
|
function "-" (Left : Complex; Right : Imaginary) return Complex is
|
|
begin
|
|
return Complex'(Left.Re, Left.Im - R (Right));
|
|
end "-";
|
|
|
|
function "-" (Left : Imaginary; Right : Complex) return Complex is
|
|
begin
|
|
return Complex'(-Right.Re, R (Left) - Right.Im);
|
|
end "-";
|
|
|
|
function "-" (Left : Imaginary; Right : Real'Base) return Complex is
|
|
begin
|
|
return Complex'(-Right, R (Left));
|
|
end "-";
|
|
|
|
function "-" (Left : Real'Base; Right : Imaginary) return Complex is
|
|
begin
|
|
return Complex'(Left, -R (Right));
|
|
end "-";
|
|
|
|
---------
|
|
-- "/" --
|
|
---------
|
|
|
|
function "/" (Left, Right : Complex) return Complex is
|
|
a : constant R := Left.Re;
|
|
b : constant R := Left.Im;
|
|
c : constant R := Right.Re;
|
|
d : constant R := Right.Im;
|
|
|
|
begin
|
|
if c = 0.0 and then d = 0.0 then
|
|
raise Constraint_Error;
|
|
else
|
|
return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
|
|
Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
|
|
end if;
|
|
end "/";
|
|
|
|
function "/" (Left, Right : Imaginary) return Real'Base is
|
|
begin
|
|
return R (Left) / R (Right);
|
|
end "/";
|
|
|
|
function "/" (Left : Complex; Right : Real'Base) return Complex is
|
|
begin
|
|
return Complex'(Left.Re / Right, Left.Im / Right);
|
|
end "/";
|
|
|
|
function "/" (Left : Real'Base; Right : Complex) return Complex is
|
|
a : constant R := Left;
|
|
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)));
|
|
end "/";
|
|
|
|
function "/" (Left : Complex; Right : Imaginary) return Complex is
|
|
a : constant R := Left.Re;
|
|
b : constant R := Left.Im;
|
|
d : constant R := R (Right);
|
|
|
|
begin
|
|
return (b / d, -(a / d));
|
|
end "/";
|
|
|
|
function "/" (Left : Imaginary; Right : Complex) return Complex is
|
|
b : constant R := R (Left);
|
|
c : constant R := Right.Re;
|
|
d : constant R := Right.Im;
|
|
|
|
begin
|
|
return (Re => b * d / (c ** 2 + d ** 2),
|
|
Im => b * c / (c ** 2 + d ** 2));
|
|
end "/";
|
|
|
|
function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
|
|
begin
|
|
return Imaginary (R (Left) / Right);
|
|
end "/";
|
|
|
|
function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
|
|
begin
|
|
return Imaginary (-(Left / R (Right)));
|
|
end "/";
|
|
|
|
---------
|
|
-- "<" --
|
|
---------
|
|
|
|
function "<" (Left, Right : Imaginary) return Boolean is
|
|
begin
|
|
return R (Left) < R (Right);
|
|
end "<";
|
|
|
|
----------
|
|
-- "<=" --
|
|
----------
|
|
|
|
function "<=" (Left, Right : Imaginary) return Boolean is
|
|
begin
|
|
return R (Left) <= R (Right);
|
|
end "<=";
|
|
|
|
---------
|
|
-- ">" --
|
|
---------
|
|
|
|
function ">" (Left, Right : Imaginary) return Boolean is
|
|
begin
|
|
return R (Left) > R (Right);
|
|
end ">";
|
|
|
|
----------
|
|
-- ">=" --
|
|
----------
|
|
|
|
function ">=" (Left, Right : Imaginary) return Boolean is
|
|
begin
|
|
return R (Left) >= R (Right);
|
|
end ">=";
|
|
|
|
-----------
|
|
-- "abs" --
|
|
-----------
|
|
|
|
function "abs" (Right : Imaginary) return Real'Base is
|
|
begin
|
|
return abs R (Right);
|
|
end "abs";
|
|
|
|
--------------
|
|
-- Argument --
|
|
--------------
|
|
|
|
function Argument (X : Complex) return Real'Base is
|
|
a : constant R := X.Re;
|
|
b : constant R := X.Im;
|
|
arg : R;
|
|
|
|
begin
|
|
if b = 0.0 then
|
|
|
|
if a >= 0.0 then
|
|
return 0.0;
|
|
else
|
|
return R'Copy_Sign (Pi, b);
|
|
end if;
|
|
|
|
elsif a = 0.0 then
|
|
|
|
if b >= 0.0 then
|
|
return Half_Pi;
|
|
else
|
|
return -Half_Pi;
|
|
end if;
|
|
|
|
else
|
|
arg := R (Atan (Double (abs (b / a))));
|
|
|
|
if a > 0.0 then
|
|
if b > 0.0 then
|
|
return arg;
|
|
else -- b < 0.0
|
|
return -arg;
|
|
end if;
|
|
|
|
else -- a < 0.0
|
|
if b >= 0.0 then
|
|
return Pi - arg;
|
|
else -- b < 0.0
|
|
return -(Pi - arg);
|
|
end if;
|
|
end if;
|
|
end if;
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
if b > 0.0 then
|
|
return Half_Pi;
|
|
else
|
|
return -Half_Pi;
|
|
end if;
|
|
end Argument;
|
|
|
|
function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
|
|
begin
|
|
if Cycle > 0.0 then
|
|
return Argument (X) * Cycle / Two_Pi;
|
|
else
|
|
raise Argument_Error;
|
|
end if;
|
|
end Argument;
|
|
|
|
----------------------------
|
|
-- Compose_From_Cartesian --
|
|
----------------------------
|
|
|
|
function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
|
|
begin
|
|
return (Re, Im);
|
|
end Compose_From_Cartesian;
|
|
|
|
function Compose_From_Cartesian (Re : Real'Base) return Complex is
|
|
begin
|
|
return (Re, 0.0);
|
|
end Compose_From_Cartesian;
|
|
|
|
function Compose_From_Cartesian (Im : Imaginary) return Complex is
|
|
begin
|
|
return (0.0, R (Im));
|
|
end Compose_From_Cartesian;
|
|
|
|
------------------------
|
|
-- Compose_From_Polar --
|
|
------------------------
|
|
|
|
function Compose_From_Polar (
|
|
Modulus, Argument : Real'Base)
|
|
return Complex
|
|
is
|
|
begin
|
|
if Modulus = 0.0 then
|
|
return (0.0, 0.0);
|
|
else
|
|
return (Modulus * R (Cos (Double (Argument))),
|
|
Modulus * R (Sin (Double (Argument))));
|
|
end if;
|
|
end Compose_From_Polar;
|
|
|
|
function Compose_From_Polar (
|
|
Modulus, Argument, Cycle : Real'Base)
|
|
return Complex
|
|
is
|
|
Arg : Real'Base;
|
|
|
|
begin
|
|
if Modulus = 0.0 then
|
|
return (0.0, 0.0);
|
|
|
|
elsif Cycle > 0.0 then
|
|
if Argument = 0.0 then
|
|
return (Modulus, 0.0);
|
|
|
|
elsif Argument = Cycle / 4.0 then
|
|
return (0.0, Modulus);
|
|
|
|
elsif Argument = Cycle / 2.0 then
|
|
return (-Modulus, 0.0);
|
|
|
|
elsif Argument = 3.0 * Cycle / R (4.0) then
|
|
return (0.0, -Modulus);
|
|
else
|
|
Arg := Two_Pi * Argument / Cycle;
|
|
return (Modulus * R (Cos (Double (Arg))),
|
|
Modulus * R (Sin (Double (Arg))));
|
|
end if;
|
|
else
|
|
raise Argument_Error;
|
|
end if;
|
|
end Compose_From_Polar;
|
|
|
|
---------------
|
|
-- Conjugate --
|
|
---------------
|
|
|
|
function Conjugate (X : Complex) return Complex is
|
|
begin
|
|
return Complex'(X.Re, -X.Im);
|
|
end Conjugate;
|
|
|
|
--------
|
|
-- Im --
|
|
--------
|
|
|
|
function Im (X : Complex) return Real'Base is
|
|
begin
|
|
return X.Im;
|
|
end Im;
|
|
|
|
function Im (X : Imaginary) return Real'Base is
|
|
begin
|
|
return R (X);
|
|
end Im;
|
|
|
|
-------------
|
|
-- Modulus --
|
|
-------------
|
|
|
|
function Modulus (X : Complex) return Real'Base is
|
|
Re2, Im2 : R;
|
|
|
|
begin
|
|
|
|
begin
|
|
Re2 := X.Re ** 2;
|
|
|
|
-- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
|
|
-- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
|
|
-- squaring does not raise constraint_error but generates infinity,
|
|
-- 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.
|
|
|
|
-- ??? same weird test, why not Re2 > R'Last ???
|
|
if not (Re2 <= R'Last) then
|
|
raise Constraint_Error;
|
|
end if;
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
return R (Double (abs (X.Re))
|
|
* Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
|
|
end;
|
|
|
|
begin
|
|
Im2 := X.Im ** 2;
|
|
|
|
-- ??? same weird test
|
|
if not (Im2 <= R'Last) then
|
|
raise Constraint_Error;
|
|
end if;
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
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
|
|
-- underflows, return the modulus of the other component. If both
|
|
-- squares underflow, use scaling as above.
|
|
|
|
if Re2 = 0.0 then
|
|
|
|
if X.Re = 0.0 then
|
|
return abs (X.Im);
|
|
|
|
elsif Im2 = 0.0 then
|
|
|
|
if X.Im = 0.0 then
|
|
return abs (X.Re);
|
|
|
|
else
|
|
if abs (X.Re) > abs (X.Im) then
|
|
return
|
|
R (Double (abs (X.Re))
|
|
* Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
|
|
else
|
|
return
|
|
R (Double (abs (X.Im))
|
|
* Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
|
|
end if;
|
|
end if;
|
|
|
|
else
|
|
return abs (X.Im);
|
|
end if;
|
|
|
|
elsif Im2 = 0.0 then
|
|
return abs (X.Re);
|
|
|
|
-- In all other cases, the naive computation will do
|
|
|
|
else
|
|
return R (Sqrt (Double (Re2 + Im2)));
|
|
end if;
|
|
end Modulus;
|
|
|
|
--------
|
|
-- Re --
|
|
--------
|
|
|
|
function Re (X : Complex) return Real'Base is
|
|
begin
|
|
return X.Re;
|
|
end Re;
|
|
|
|
------------
|
|
-- Set_Im --
|
|
------------
|
|
|
|
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 : Real'Base) is
|
|
begin
|
|
X := Imaginary (Im);
|
|
end Set_Im;
|
|
|
|
------------
|
|
-- Set_Re --
|
|
------------
|
|
|
|
procedure Set_Re (X : in out Complex; Re : Real'Base) is
|
|
begin
|
|
X.Re := Re;
|
|
end Set_Re;
|
|
|
|
end Ada.Numerics.Generic_Complex_Types;
|