diff options
Diffstat (limited to 'gcc-4.8.3/gcc/ada/a-ngcefu.adb')
| -rw-r--r-- | gcc-4.8.3/gcc/ada/a-ngcefu.adb | 708 |
1 files changed, 708 insertions, 0 deletions
diff --git a/gcc-4.8.3/gcc/ada/a-ngcefu.adb b/gcc-4.8.3/gcc/ada/a-ngcefu.adb new file mode 100644 index 000000000..edcdb5a72 --- /dev/null +++ b/gcc-4.8.3/gcc/ada/a-ngcefu.adb @@ -0,0 +1,708 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 1992-2009, 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; |