711 lines
18 KiB
Ada
711 lines
18 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT RUN-TIME COMPONENTS --
|
|
-- --
|
|
-- ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 1992-2013, 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.Generic_Elementary_Functions;
|
|
|
|
package body Ada.Numerics.Generic_Complex_Elementary_Functions is
|
|
|
|
package Elementary_Functions is new
|
|
Ada.Numerics.Generic_Elementary_Functions (Real'Base);
|
|
use Elementary_Functions;
|
|
|
|
PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
|
|
PI_2 : constant := PI / 2.0;
|
|
Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
|
|
Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
|
|
|
|
subtype T is Real'Base;
|
|
|
|
Epsilon : constant T := 2.0 ** (1 - T'Model_Mantissa);
|
|
Square_Root_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
|
|
Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
|
|
Root_Root_Epsilon : constant T := Sqrt_Two **
|
|
((1 - T'Model_Mantissa) / 2);
|
|
Log_Inverse_Epsilon_2 : constant T := T (T'Model_Mantissa - 1) / 2.0;
|
|
|
|
Complex_Zero : constant Complex := (0.0, 0.0);
|
|
Complex_One : constant Complex := (1.0, 0.0);
|
|
Complex_I : constant Complex := (0.0, 1.0);
|
|
Half_Pi : constant Complex := (PI_2, 0.0);
|
|
|
|
--------
|
|
-- ** --
|
|
--------
|
|
|
|
function "**" (Left : Complex; Right : Complex) return Complex is
|
|
begin
|
|
if Re (Right) = 0.0
|
|
and then Im (Right) = 0.0
|
|
and then Re (Left) = 0.0
|
|
and then Im (Left) = 0.0
|
|
then
|
|
raise Argument_Error;
|
|
|
|
elsif Re (Left) = 0.0
|
|
and then Im (Left) = 0.0
|
|
and then Re (Right) < 0.0
|
|
then
|
|
raise Constraint_Error;
|
|
|
|
elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
|
|
return Left;
|
|
|
|
elsif Right = (0.0, 0.0) then
|
|
return Complex_One;
|
|
|
|
elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
|
|
return 1.0 + Right;
|
|
|
|
elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
|
|
return Left;
|
|
|
|
else
|
|
return Exp (Right * Log (Left));
|
|
end if;
|
|
end "**";
|
|
|
|
function "**" (Left : Real'Base; Right : Complex) return Complex is
|
|
begin
|
|
if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
|
|
raise Argument_Error;
|
|
|
|
elsif Left = 0.0 and then Re (Right) < 0.0 then
|
|
raise Constraint_Error;
|
|
|
|
elsif Left = 0.0 then
|
|
return Compose_From_Cartesian (Left, 0.0);
|
|
|
|
elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
|
|
return Complex_One;
|
|
|
|
elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
|
|
return Compose_From_Cartesian (Left, 0.0);
|
|
|
|
else
|
|
return Exp (Log (Left) * Right);
|
|
end if;
|
|
end "**";
|
|
|
|
function "**" (Left : Complex; Right : Real'Base) return Complex is
|
|
begin
|
|
if Right = 0.0
|
|
and then Re (Left) = 0.0
|
|
and then Im (Left) = 0.0
|
|
then
|
|
raise Argument_Error;
|
|
|
|
elsif Re (Left) = 0.0
|
|
and then Im (Left) = 0.0
|
|
and then Right < 0.0
|
|
then
|
|
raise Constraint_Error;
|
|
|
|
elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
|
|
return Left;
|
|
|
|
elsif Right = 0.0 then
|
|
return Complex_One;
|
|
|
|
elsif Right = 1.0 then
|
|
return Left;
|
|
|
|
else
|
|
return Exp (Right * Log (Left));
|
|
end if;
|
|
end "**";
|
|
|
|
------------
|
|
-- Arccos --
|
|
------------
|
|
|
|
function Arccos (X : Complex) return Complex is
|
|
Result : Complex;
|
|
|
|
begin
|
|
if X = Complex_One then
|
|
return Complex_Zero;
|
|
|
|
elsif abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return Half_Pi - X;
|
|
|
|
elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
abs Im (X) > Inv_Square_Root_Epsilon
|
|
then
|
|
return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
|
|
Complex_I * Sqrt ((1.0 - X) / 2.0));
|
|
end if;
|
|
|
|
Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
|
|
|
|
if Im (X) = 0.0
|
|
and then abs Re (X) <= 1.00
|
|
then
|
|
Set_Im (Result, Im (X));
|
|
end if;
|
|
|
|
return Result;
|
|
end Arccos;
|
|
|
|
-------------
|
|
-- Arccosh --
|
|
-------------
|
|
|
|
function Arccosh (X : Complex) return Complex is
|
|
Result : Complex;
|
|
|
|
begin
|
|
if X = Complex_One then
|
|
return Complex_Zero;
|
|
|
|
elsif abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
|
|
|
|
elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
abs Im (X) > Inv_Square_Root_Epsilon
|
|
then
|
|
Result := Log_Two + Log (X);
|
|
|
|
else
|
|
Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
|
|
Sqrt ((X - 1.0) / 2.0));
|
|
end if;
|
|
|
|
if Re (Result) <= 0.0 then
|
|
Result := -Result;
|
|
end if;
|
|
|
|
return Result;
|
|
end Arccosh;
|
|
|
|
------------
|
|
-- Arccot --
|
|
------------
|
|
|
|
function Arccot (X : Complex) return Complex is
|
|
Xt : Complex;
|
|
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return Half_Pi - X;
|
|
|
|
elsif abs Re (X) > 1.0 / Epsilon or else
|
|
abs Im (X) > 1.0 / Epsilon
|
|
then
|
|
Xt := Complex_One / X;
|
|
|
|
if Re (X) < 0.0 then
|
|
Set_Re (Xt, PI - Re (Xt));
|
|
return Xt;
|
|
else
|
|
return Xt;
|
|
end if;
|
|
end if;
|
|
|
|
Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
|
|
|
|
if Re (Xt) < 0.0 then
|
|
Xt := PI + Xt;
|
|
end if;
|
|
|
|
return Xt;
|
|
end Arccot;
|
|
|
|
--------------
|
|
-- Arccoth --
|
|
--------------
|
|
|
|
function Arccoth (X : Complex) return Complex is
|
|
R : Complex;
|
|
|
|
begin
|
|
if X = (0.0, 0.0) then
|
|
return Compose_From_Cartesian (0.0, PI_2);
|
|
|
|
elsif abs Re (X) < Square_Root_Epsilon
|
|
and then abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return PI_2 * Complex_I + X;
|
|
|
|
elsif abs Re (X) > 1.0 / Epsilon or else
|
|
abs Im (X) > 1.0 / Epsilon
|
|
then
|
|
if Im (X) > 0.0 then
|
|
return (0.0, 0.0);
|
|
else
|
|
return PI * Complex_I;
|
|
end if;
|
|
|
|
elsif Im (X) = 0.0 and then Re (X) = 1.0 then
|
|
raise Constraint_Error;
|
|
|
|
elsif Im (X) = 0.0 and then Re (X) = -1.0 then
|
|
raise Constraint_Error;
|
|
end if;
|
|
|
|
begin
|
|
R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
|
|
end;
|
|
|
|
if Im (R) < 0.0 then
|
|
Set_Im (R, PI + Im (R));
|
|
end if;
|
|
|
|
if Re (X) = 0.0 then
|
|
Set_Re (R, Re (X));
|
|
end if;
|
|
|
|
return R;
|
|
end Arccoth;
|
|
|
|
------------
|
|
-- Arcsin --
|
|
------------
|
|
|
|
function Arcsin (X : Complex) return Complex is
|
|
Result : Complex;
|
|
|
|
begin
|
|
-- For very small argument, sin (x) = x
|
|
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
|
|
elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
abs Im (X) > Inv_Square_Root_Epsilon
|
|
then
|
|
Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
|
|
|
|
if Im (Result) > PI_2 then
|
|
Set_Im (Result, PI - Im (X));
|
|
|
|
elsif Im (Result) < -PI_2 then
|
|
Set_Im (Result, -(PI + Im (X)));
|
|
end if;
|
|
|
|
return Result;
|
|
end if;
|
|
|
|
Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
|
|
|
|
if Re (X) = 0.0 then
|
|
Set_Re (Result, Re (X));
|
|
|
|
elsif Im (X) = 0.0
|
|
and then abs Re (X) <= 1.00
|
|
then
|
|
Set_Im (Result, Im (X));
|
|
end if;
|
|
|
|
return Result;
|
|
end Arcsin;
|
|
|
|
-------------
|
|
-- Arcsinh --
|
|
-------------
|
|
|
|
function Arcsinh (X : Complex) return Complex is
|
|
Result : Complex;
|
|
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
|
|
elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
abs Im (X) > Inv_Square_Root_Epsilon
|
|
then
|
|
Result := Log_Two + Log (X); -- may have wrong sign
|
|
|
|
if (Re (X) < 0.0 and then Re (Result) > 0.0)
|
|
or else (Re (X) > 0.0 and then Re (Result) < 0.0)
|
|
then
|
|
Set_Re (Result, -Re (Result));
|
|
end if;
|
|
|
|
return Result;
|
|
end if;
|
|
|
|
Result := Log (X + Sqrt (1.0 + X * X));
|
|
|
|
if Re (X) = 0.0 then
|
|
Set_Re (Result, Re (X));
|
|
elsif Im (X) = 0.0 then
|
|
Set_Im (Result, Im (X));
|
|
end if;
|
|
|
|
return Result;
|
|
end Arcsinh;
|
|
|
|
------------
|
|
-- Arctan --
|
|
------------
|
|
|
|
function Arctan (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
|
|
else
|
|
return -Complex_I * (Log (1.0 + Complex_I * X)
|
|
- Log (1.0 - Complex_I * X)) / 2.0;
|
|
end if;
|
|
end Arctan;
|
|
|
|
-------------
|
|
-- Arctanh --
|
|
-------------
|
|
|
|
function Arctanh (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
else
|
|
return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
|
|
end if;
|
|
end Arctanh;
|
|
|
|
---------
|
|
-- Cos --
|
|
---------
|
|
|
|
function Cos (X : Complex) return Complex is
|
|
begin
|
|
return
|
|
Compose_From_Cartesian
|
|
(Cos (Re (X)) * Cosh (Im (X)),
|
|
-(Sin (Re (X)) * Sinh (Im (X))));
|
|
end Cos;
|
|
|
|
----------
|
|
-- Cosh --
|
|
----------
|
|
|
|
function Cosh (X : Complex) return Complex is
|
|
begin
|
|
return
|
|
Compose_From_Cartesian
|
|
(Cosh (Re (X)) * Cos (Im (X)),
|
|
Sinh (Re (X)) * Sin (Im (X)));
|
|
end Cosh;
|
|
|
|
---------
|
|
-- Cot --
|
|
---------
|
|
|
|
function Cot (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return Complex_One / X;
|
|
|
|
elsif Im (X) > Log_Inverse_Epsilon_2 then
|
|
return -Complex_I;
|
|
|
|
elsif Im (X) < -Log_Inverse_Epsilon_2 then
|
|
return Complex_I;
|
|
end if;
|
|
|
|
return Cos (X) / Sin (X);
|
|
end Cot;
|
|
|
|
----------
|
|
-- Coth --
|
|
----------
|
|
|
|
function Coth (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return Complex_One / X;
|
|
|
|
elsif Re (X) > Log_Inverse_Epsilon_2 then
|
|
return Complex_One;
|
|
|
|
elsif Re (X) < -Log_Inverse_Epsilon_2 then
|
|
return -Complex_One;
|
|
|
|
else
|
|
return Cosh (X) / Sinh (X);
|
|
end if;
|
|
end Coth;
|
|
|
|
---------
|
|
-- Exp --
|
|
---------
|
|
|
|
function Exp (X : Complex) return Complex is
|
|
EXP_RE_X : constant Real'Base := Exp (Re (X));
|
|
|
|
begin
|
|
return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
|
|
EXP_RE_X * Sin (Im (X)));
|
|
end Exp;
|
|
|
|
function Exp (X : Imaginary) return Complex is
|
|
ImX : constant Real'Base := Im (X);
|
|
|
|
begin
|
|
return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
|
|
end Exp;
|
|
|
|
---------
|
|
-- Log --
|
|
---------
|
|
|
|
function Log (X : Complex) return Complex is
|
|
ReX : Real'Base;
|
|
ImX : Real'Base;
|
|
Z : Complex;
|
|
|
|
begin
|
|
if Re (X) = 0.0 and then Im (X) = 0.0 then
|
|
raise Constraint_Error;
|
|
|
|
elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
|
|
and then abs Im (X) < Root_Root_Epsilon
|
|
then
|
|
Z := X;
|
|
Set_Re (Z, Re (Z) - 1.0);
|
|
|
|
return (1.0 - (1.0 / 2.0 -
|
|
(1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
|
|
end if;
|
|
|
|
begin
|
|
ReX := Log (Modulus (X));
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
ReX := Log (Modulus (X / 2.0)) - Log_Two;
|
|
end;
|
|
|
|
ImX := Arctan (Im (X), Re (X));
|
|
|
|
if ImX > PI then
|
|
ImX := ImX - 2.0 * PI;
|
|
end if;
|
|
|
|
return Compose_From_Cartesian (ReX, ImX);
|
|
end Log;
|
|
|
|
---------
|
|
-- Sin --
|
|
---------
|
|
|
|
function Sin (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon
|
|
and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
end if;
|
|
|
|
return
|
|
Compose_From_Cartesian
|
|
(Sin (Re (X)) * Cosh (Im (X)),
|
|
Cos (Re (X)) * Sinh (Im (X)));
|
|
end Sin;
|
|
|
|
----------
|
|
-- Sinh --
|
|
----------
|
|
|
|
function Sinh (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
|
|
else
|
|
return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
|
|
Cosh (Re (X)) * Sin (Im (X)));
|
|
end if;
|
|
end Sinh;
|
|
|
|
----------
|
|
-- Sqrt --
|
|
----------
|
|
|
|
function Sqrt (X : Complex) return Complex is
|
|
ReX : constant Real'Base := Re (X);
|
|
ImX : constant Real'Base := Im (X);
|
|
XR : constant Real'Base := abs Re (X);
|
|
YR : constant Real'Base := abs Im (X);
|
|
R : Real'Base;
|
|
R_X : Real'Base;
|
|
R_Y : Real'Base;
|
|
|
|
begin
|
|
-- Deal with pure real case, see (RM G.1.2(39))
|
|
|
|
if ImX = 0.0 then
|
|
if ReX > 0.0 then
|
|
return
|
|
Compose_From_Cartesian
|
|
(Sqrt (ReX), 0.0);
|
|
|
|
elsif ReX = 0.0 then
|
|
return X;
|
|
|
|
else
|
|
return
|
|
Compose_From_Cartesian
|
|
(0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
|
|
end if;
|
|
|
|
elsif ReX = 0.0 then
|
|
R_X := Sqrt (YR / 2.0);
|
|
|
|
if ImX > 0.0 then
|
|
return Compose_From_Cartesian (R_X, R_X);
|
|
else
|
|
return Compose_From_Cartesian (R_X, -R_X);
|
|
end if;
|
|
|
|
else
|
|
R := Sqrt (XR ** 2 + YR ** 2);
|
|
|
|
-- If the square of the modulus overflows, try rescaling the
|
|
-- real and imaginary parts. We cannot depend on an exception
|
|
-- being raised on all targets.
|
|
|
|
if R > Real'Base'Last then
|
|
raise Constraint_Error;
|
|
end if;
|
|
|
|
-- We are solving the system
|
|
|
|
-- XR = R_X ** 2 - Y_R ** 2 (1)
|
|
-- YR = 2.0 * R_X * R_Y (2)
|
|
--
|
|
-- The symmetric solution involves square roots for both R_X and
|
|
-- R_Y, but it is more accurate to use the square root with the
|
|
-- larger argument for either R_X or R_Y, and equation (2) for the
|
|
-- other.
|
|
|
|
if ReX < 0.0 then
|
|
R_Y := Sqrt (0.5 * (R - ReX));
|
|
R_X := YR / (2.0 * R_Y);
|
|
|
|
else
|
|
R_X := Sqrt (0.5 * (R + ReX));
|
|
R_Y := YR / (2.0 * R_X);
|
|
end if;
|
|
end if;
|
|
|
|
if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
|
|
R_Y := -R_Y;
|
|
end if;
|
|
return Compose_From_Cartesian (R_X, R_Y);
|
|
|
|
exception
|
|
when Constraint_Error =>
|
|
|
|
-- Rescale and try again
|
|
|
|
R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
|
|
R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
|
|
R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
|
|
|
|
if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
|
|
R_Y := -R_Y;
|
|
end if;
|
|
|
|
return Compose_From_Cartesian (R_X, R_Y);
|
|
end Sqrt;
|
|
|
|
---------
|
|
-- Tan --
|
|
---------
|
|
|
|
function Tan (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
|
|
elsif Im (X) > Log_Inverse_Epsilon_2 then
|
|
return Complex_I;
|
|
|
|
elsif Im (X) < -Log_Inverse_Epsilon_2 then
|
|
return -Complex_I;
|
|
|
|
else
|
|
return Sin (X) / Cos (X);
|
|
end if;
|
|
end Tan;
|
|
|
|
----------
|
|
-- Tanh --
|
|
----------
|
|
|
|
function Tanh (X : Complex) return Complex is
|
|
begin
|
|
if abs Re (X) < Square_Root_Epsilon and then
|
|
abs Im (X) < Square_Root_Epsilon
|
|
then
|
|
return X;
|
|
|
|
elsif Re (X) > Log_Inverse_Epsilon_2 then
|
|
return Complex_One;
|
|
|
|
elsif Re (X) < -Log_Inverse_Epsilon_2 then
|
|
return -Complex_One;
|
|
|
|
else
|
|
return Sinh (X) / Cosh (X);
|
|
end if;
|
|
end Tanh;
|
|
|
|
end Ada.Numerics.Generic_Complex_Elementary_Functions;
|