diff options
Diffstat (limited to 'gcc-4.4.0/gcc/ada/math_lib.adb')
| -rw-r--r-- | gcc-4.4.0/gcc/ada/math_lib.adb | 1025 |
1 files changed, 0 insertions, 1025 deletions
diff --git a/gcc-4.4.0/gcc/ada/math_lib.adb b/gcc-4.4.0/gcc/ada/math_lib.adb deleted file mode 100644 index e539f477b..000000000 --- a/gcc-4.4.0/gcc/ada/math_lib.adb +++ /dev/null @@ -1,1025 +0,0 @@ ------------------------------------------------------------------------------- --- -- --- GNAT RUN-TIME COMPONENTS -- --- -- --- M A T H _ L I B -- --- -- --- 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. -- --- -- ------------------------------------------------------------------------------- - --- This body is specifically for using an Ada interface to C math.h to get --- the computation engine. Many special cases are handled locally to avoid --- unnecessary calls. This is not a "strict" implementation, but takes full --- advantage of the C functions, e.g. in providing interface to hardware --- provided versions of the elementary functions. - --- A known weakness is that on the x86, all computation is done in Double, --- which means that a lot of accuracy is lost for the Long_Long_Float case. - --- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, --- sinh, cosh, tanh from C library via math.h - --- This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that --- provides a compatible body for the DEC Math_Lib package. - -with Ada.Numerics.Aux; -use type Ada.Numerics.Aux.Double; -with Ada.Numerics; use Ada.Numerics; - -package body Math_Lib is - - Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755; - - Two_Pi : constant Real'Base := 2.0 * Pi; - Half_Pi : constant Real'Base := Pi / 2.0; - Fourth_Pi : constant Real'Base := Pi / 4.0; - Epsilon : constant Real'Base := Real'Base'Epsilon; - IEpsilon : constant Real'Base := 1.0 / Epsilon; - - subtype Double is Aux.Double; - - DEpsilon : constant Double := Double (Epsilon); - DIEpsilon : constant Double := Double (IEpsilon); - - ----------------------- - -- Local Subprograms -- - ----------------------- - - function Arctan - (Y : Real; - A : Real := 1.0) - return Real; - - function Arctan - (Y : Real; - A : Real := 1.0; - Cycle : Real) - return Real; - - function Exact_Remainder - (A : Real; - Y : Real) - return Real; - -- Computes exact remainder of A divided by Y - - function Half_Log_Epsilon return Real; - -- Function to provide constant: 0.5 * Log (Epsilon) - - function Local_Atan - (Y : Real; - A : Real := 1.0) - return Real; - -- Common code for arc tangent after cycle reduction - - function Log_Inverse_Epsilon return Real; - -- Function to provide constant: Log (1.0 / Epsilon) - - function Square_Root_Epsilon return Real; - -- Function to provide constant: Sqrt (Epsilon) - - ---------- - -- "**" -- - ---------- - - function "**" (A1, A2 : Real) return Real is - - begin - if A1 = 0.0 - and then A2 = 0.0 - then - raise Argument_Error; - - elsif A1 < 0.0 then - raise Argument_Error; - - elsif A2 = 0.0 then - return 1.0; - - elsif A1 = 0.0 then - if A2 < 0.0 then - raise Constraint_Error; - else - return 0.0; - end if; - - elsif A1 = 1.0 then - return 1.0; - - elsif A2 = 1.0 then - return A1; - - else - begin - if A2 = 2.0 then - return A1 * A1; - else - return - Real (Aux.pow (Double (A1), Double (A2))); - end if; - - exception - when others => - raise Constraint_Error; - end; - end if; - end "**"; - - ------------ - -- Arccos -- - ------------ - - -- Natural cycle - - function Arccos (A : Real) return Real is - Temp : Real'Base; - - begin - if abs A > 1.0 then - raise Argument_Error; - - elsif abs A < Square_Root_Epsilon then - return Pi / 2.0 - A; - - elsif A = 1.0 then - return 0.0; - - elsif A = -1.0 then - return Pi; - end if; - - Temp := Real (Aux.acos (Double (A))); - - if Temp < 0.0 then - Temp := Pi + Temp; - end if; - - return Temp; - end Arccos; - - -- Arbitrary cycle - - function Arccos (A, Cycle : Real) return Real is - Temp : Real'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif abs A > 1.0 then - raise Argument_Error; - - elsif abs A < Square_Root_Epsilon then - return Cycle / 4.0; - - elsif A = 1.0 then - return 0.0; - - elsif A = -1.0 then - return Cycle / 2.0; - end if; - - Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle); - - if Temp < 0.0 then - Temp := Cycle / 2.0 + Temp; - end if; - - return Temp; - end Arccos; - - ------------- - -- Arccosh -- - ------------- - - function Arccosh (A : Real) return Real is - begin - -- Return Log (A - Sqrt (A * A - 1.0)); double valued, - -- only positive value returned - -- What is this comment ??? - - if A < 1.0 then - raise Argument_Error; - - elsif A < 1.0 + Square_Root_Epsilon then - return A - 1.0; - - elsif abs A > 1.0 / Square_Root_Epsilon then - return Log (A) + Log_Two; - - else - return Log (A + Sqrt (A * A - 1.0)); - end if; - end Arccosh; - - ------------ - -- Arccot -- - ------------ - - -- Natural cycle - - function Arccot - (A : Real; - Y : Real := 1.0) - return Real - is - begin - -- Just reverse arguments - - return Arctan (Y, A); - end Arccot; - - -- Arbitrary cycle - - function Arccot - (A : Real; - Y : Real := 1.0; - Cycle : Real) - return Real - is - begin - -- Just reverse arguments - - return Arctan (Y, A, Cycle); - end Arccot; - - ------------- - -- Arccoth -- - ------------- - - function Arccoth (A : Real) return Real is - begin - if abs A = 1.0 then - raise Constraint_Error; - - elsif abs A < 1.0 then - raise Argument_Error; - - elsif abs A > 1.0 / Epsilon then - return 0.0; - - else - return 0.5 * Log ((1.0 + A) / (A - 1.0)); - end if; - end Arccoth; - - ------------ - -- Arcsin -- - ------------ - - -- Natural cycle - - function Arcsin (A : Real) return Real is - begin - if abs A > 1.0 then - raise Argument_Error; - - elsif abs A < Square_Root_Epsilon then - return A; - - elsif A = 1.0 then - return Pi / 2.0; - - elsif A = -1.0 then - return -Pi / 2.0; - end if; - - return Real (Aux.asin (Double (A))); - end Arcsin; - - -- Arbitrary cycle - - function Arcsin (A, Cycle : Real) return Real is - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif abs A > 1.0 then - raise Argument_Error; - - elsif A = 0.0 then - return A; - - elsif A = 1.0 then - return Cycle / 4.0; - - elsif A = -1.0 then - return -Cycle / 4.0; - end if; - - return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle); - end Arcsin; - - ------------- - -- Arcsinh -- - ------------- - - function Arcsinh (A : Real) return Real is - begin - if abs A < Square_Root_Epsilon then - return A; - - elsif A > 1.0 / Square_Root_Epsilon then - return Log (A) + Log_Two; - - elsif A < -1.0 / Square_Root_Epsilon then - return -(Log (-A) + Log_Two); - - elsif A < 0.0 then - return -Log (abs A + Sqrt (A * A + 1.0)); - - else - return Log (A + Sqrt (A * A + 1.0)); - end if; - end Arcsinh; - - ------------ - -- Arctan -- - ------------ - - -- Natural cycle - - function Arctan - (Y : Real; - A : Real := 1.0) - return Real - is - begin - if A = 0.0 - and then Y = 0.0 - then - raise Argument_Error; - - elsif Y = 0.0 then - if A > 0.0 then - return 0.0; - else -- A < 0.0 - return Pi; - end if; - - elsif A = 0.0 then - if Y > 0.0 then - return Half_Pi; - else -- Y < 0.0 - return -Half_Pi; - end if; - - else - return Local_Atan (Y, A); - end if; - end Arctan; - - -- Arbitrary cycle - - function Arctan - (Y : Real; - A : Real := 1.0; - Cycle : Real) - return Real - is - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif A = 0.0 - and then Y = 0.0 - then - raise Argument_Error; - - elsif Y = 0.0 then - if A > 0.0 then - return 0.0; - else -- A < 0.0 - return Cycle / 2.0; - end if; - - elsif A = 0.0 then - if Y > 0.0 then - return Cycle / 4.0; - else -- Y < 0.0 - return -Cycle / 4.0; - end if; - - else - return Local_Atan (Y, A) * Cycle / Two_Pi; - end if; - end Arctan; - - ------------- - -- Arctanh -- - ------------- - - function Arctanh (A : Real) return Real is - begin - if abs A = 1.0 then - raise Constraint_Error; - - elsif abs A > 1.0 then - raise Argument_Error; - - elsif abs A < Square_Root_Epsilon then - return A; - - else - return 0.5 * Log ((1.0 + A) / (1.0 - A)); - end if; - end Arctanh; - - --------- - -- Cos -- - --------- - - -- Natural cycle - - function Cos (A : Real) return Real is - begin - if A = 0.0 then - return 1.0; - - elsif abs A < Square_Root_Epsilon then - return 1.0; - - end if; - - return Real (Aux.Cos (Double (A))); - end Cos; - - -- Arbitrary cycle - - function Cos (A, Cycle : Real) return Real is - T : Real'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif A = 0.0 then - return 1.0; - end if; - - T := Exact_Remainder (abs (A), Cycle) / Cycle; - - if T = 0.25 - or else T = 0.75 - or else T = -0.25 - or else T = -0.75 - then - return 0.0; - - elsif T = 0.5 or T = -0.5 then - return -1.0; - end if; - - return Real (Aux.Cos (Double (T * Two_Pi))); - end Cos; - - ---------- - -- Cosh -- - ---------- - - function Cosh (A : Real) return Real is - begin - if abs A < Square_Root_Epsilon then - return 1.0; - - elsif abs A > Log_Inverse_Epsilon then - return Exp ((abs A) - Log_Two); - end if; - - return Real (Aux.cosh (Double (A))); - - exception - when others => - raise Constraint_Error; - end Cosh; - - --------- - -- Cot -- - --------- - - -- Natural cycle - - function Cot (A : Real) return Real is - begin - if A = 0.0 then - raise Constraint_Error; - - elsif abs A < Square_Root_Epsilon then - return 1.0 / A; - end if; - - return Real (1.0 / Real'Base (Aux.tan (Double (A)))); - end Cot; - - -- Arbitrary cycle - - function Cot (A, Cycle : Real) return Real is - T : Real'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif A = 0.0 then - raise Constraint_Error; - - elsif abs A < Square_Root_Epsilon then - return 1.0 / A; - end if; - - T := Exact_Remainder (A, Cycle) / Cycle; - - if T = 0.0 or T = 0.5 or T = -0.5 then - raise Constraint_Error; - else - return Cos (T * Two_Pi) / Sin (T * Two_Pi); - end if; - end Cot; - - ---------- - -- Coth -- - ---------- - - function Coth (A : Real) return Real is - begin - if A = 0.0 then - raise Constraint_Error; - - elsif A < Half_Log_Epsilon then - return -1.0; - - elsif A > -Half_Log_Epsilon then - return 1.0; - - elsif abs A < Square_Root_Epsilon then - return 1.0 / A; - end if; - - return Real (1.0 / Real'Base (Aux.tanh (Double (A)))); - end Coth; - - --------------------- - -- Exact_Remainder -- - --------------------- - - function Exact_Remainder - (A : Real; - Y : Real) - return Real - is - Denominator : Real'Base := abs A; - Divisor : Real'Base := abs Y; - Reducer : Real'Base; - Sign : Real'Base := 1.0; - - begin - if Y = 0.0 then - raise Constraint_Error; - - elsif A = 0.0 then - return 0.0; - - elsif A = Y then - return 0.0; - - elsif Denominator < Divisor then - return A; - end if; - - while Denominator >= Divisor loop - - -- Put divisors mantissa with denominators exponent to make reducer - - Reducer := Divisor; - - begin - while Reducer * 1_048_576.0 < Denominator loop - Reducer := Reducer * 1_048_576.0; - end loop; - - exception - when others => null; - end; - - begin - while Reducer * 1_024.0 < Denominator loop - Reducer := Reducer * 1_024.0; - end loop; - - exception - when others => null; - end; - - begin - while Reducer * 2.0 < Denominator loop - Reducer := Reducer * 2.0; - end loop; - - exception - when others => null; - end; - - Denominator := Denominator - Reducer; - end loop; - - if A < 0.0 then - return -Denominator; - else - return Denominator; - end if; - end Exact_Remainder; - - --------- - -- Exp -- - --------- - - function Exp (A : Real) return Real is - Result : Real'Base; - - begin - if A = 0.0 then - return 1.0; - - else - Result := Real (Aux.Exp (Double (A))); - - -- The check here catches the case of Exp returning IEEE infinity - - if Result > Real'Last then - raise Constraint_Error; - else - return Result; - end if; - end if; - end Exp; - - ---------------------- - -- Half_Log_Epsilon -- - ---------------------- - - -- Cannot precompute this constant, because this is required to be a - -- pure package, which allows no state. A pity, but no way around it! - - function Half_Log_Epsilon return Real is - begin - return Real (0.5 * Real'Base (Aux.Log (DEpsilon))); - end Half_Log_Epsilon; - - ---------------- - -- Local_Atan -- - ---------------- - - function Local_Atan - (Y : Real; - A : Real := 1.0) - return Real - is - Z : Real'Base; - Raw_Atan : Real'Base; - - begin - if abs Y > abs A then - Z := abs (A / Y); - else - Z := abs (Y / A); - end if; - - if Z < Square_Root_Epsilon then - Raw_Atan := Z; - - elsif Z = 1.0 then - Raw_Atan := Pi / 4.0; - - elsif Z < Square_Root_Epsilon then - Raw_Atan := Z; - - else - Raw_Atan := Real'Base (Aux.Atan (Double (Z))); - end if; - - if abs Y > abs A then - Raw_Atan := Half_Pi - Raw_Atan; - end if; - - if A > 0.0 then - if Y > 0.0 then - return Raw_Atan; - else -- Y < 0.0 - return -Raw_Atan; - end if; - - else -- A < 0.0 - if Y > 0.0 then - return Pi - Raw_Atan; - else -- Y < 0.0 - return -(Pi - Raw_Atan); - end if; - end if; - end Local_Atan; - - --------- - -- Log -- - --------- - - -- Natural base - - function Log (A : Real) return Real is - begin - if A < 0.0 then - raise Argument_Error; - - elsif A = 0.0 then - raise Constraint_Error; - - elsif A = 1.0 then - return 0.0; - end if; - - return Real (Aux.Log (Double (A))); - end Log; - - -- Arbitrary base - - function Log (A, Base : Real) return Real is - begin - if A < 0.0 then - raise Argument_Error; - - elsif Base <= 0.0 or else Base = 1.0 then - raise Argument_Error; - - elsif A = 0.0 then - raise Constraint_Error; - - elsif A = 1.0 then - return 0.0; - end if; - - return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base))); - end Log; - - ------------------------- - -- Log_Inverse_Epsilon -- - ------------------------- - - -- Cannot precompute this constant, because this is required to be a - -- pure package, which allows no state. A pity, but no way around it! - - function Log_Inverse_Epsilon return Real is - begin - return Real (Aux.Log (DIEpsilon)); - end Log_Inverse_Epsilon; - - --------- - -- Sin -- - --------- - - -- Natural cycle - - function Sin (A : Real) return Real is - begin - if abs A < Square_Root_Epsilon then - return A; - end if; - - return Real (Aux.Sin (Double (A))); - end Sin; - - -- Arbitrary cycle - - function Sin (A, Cycle : Real) return Real is - T : Real'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif A = 0.0 then - return A; - end if; - - T := Exact_Remainder (A, Cycle) / Cycle; - - if T = 0.0 or T = 0.5 or T = -0.5 then - return 0.0; - - elsif T = 0.25 or T = -0.75 then - return 1.0; - - elsif T = -0.25 or T = 0.75 then - return -1.0; - - end if; - - return Real (Aux.Sin (Double (T * Two_Pi))); - end Sin; - - ---------- - -- Sinh -- - ---------- - - function Sinh (A : Real) return Real is - begin - if abs A < Square_Root_Epsilon then - return A; - - elsif A > Log_Inverse_Epsilon then - return Exp (A - Log_Two); - - elsif A < -Log_Inverse_Epsilon then - return -Exp ((-A) - Log_Two); - end if; - - return Real (Aux.Sinh (Double (A))); - - exception - when others => - raise Constraint_Error; - end Sinh; - - ------------------------- - -- Square_Root_Epsilon -- - ------------------------- - - -- Cannot precompute this constant, because this is required to be a - -- pure package, which allows no state. A pity, but no way around it! - - function Square_Root_Epsilon return Real is - begin - return Real (Aux.Sqrt (DEpsilon)); - end Square_Root_Epsilon; - - ---------- - -- Sqrt -- - ---------- - - function Sqrt (A : Real) return Real is - begin - if A < 0.0 then - raise Argument_Error; - - -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE - - elsif A = 0.0 then - return A; - - -- Sqrt (1.0) must be exact for good complex accuracy - - elsif A = 1.0 then - return 1.0; - - end if; - - return Real (Aux.Sqrt (Double (A))); - end Sqrt; - - --------- - -- Tan -- - --------- - - -- Natural cycle - - function Tan (A : Real) return Real is - begin - if abs A < Square_Root_Epsilon then - return A; - - elsif abs A = Pi / 2.0 then - raise Constraint_Error; - end if; - - return Real (Aux.tan (Double (A))); - end Tan; - - -- Arbitrary cycle - - function Tan (A, Cycle : Real) return Real is - T : Real'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif A = 0.0 then - return A; - end if; - - T := Exact_Remainder (A, Cycle) / Cycle; - - if T = 0.25 - or else T = 0.75 - or else T = -0.25 - or else T = -0.75 - then - raise Constraint_Error; - - else - return Sin (T * Two_Pi) / Cos (T * Two_Pi); - end if; - end Tan; - - ---------- - -- Tanh -- - ---------- - - function Tanh (A : Real) return Real is - begin - if A < Half_Log_Epsilon then - return -1.0; - - elsif A > -Half_Log_Epsilon then - return 1.0; - - elsif abs A < Square_Root_Epsilon then - return A; - end if; - - return Real (Aux.tanh (Double (A))); - end Tanh; - - ---------------------------- - -- DEC-Specific functions -- - ---------------------------- - - function LOG10 (A : REAL) return REAL is - begin - return Log (A, 10.0); - end LOG10; - - function LOG2 (A : REAL) return REAL is - begin - return Log (A, 2.0); - end LOG2; - - function ASIN (A : REAL) return REAL renames Arcsin; - function ACOS (A : REAL) return REAL renames Arccos; - - function ATAN (A : REAL) return REAL is - begin - return Arctan (A, 1.0); - end ATAN; - - function ATAN2 (A1, A2 : REAL) return REAL renames Arctan; - - function SIND (A : REAL) return REAL is - begin - return Sin (A, 360.0); - end SIND; - - function COSD (A : REAL) return REAL is - begin - return Cos (A, 360.0); - end COSD; - - function TAND (A : REAL) return REAL is - begin - return Tan (A, 360.0); - end TAND; - - function ASIND (A : REAL) return REAL is - begin - return Arcsin (A, 360.0); - end ASIND; - - function ACOSD (A : REAL) return REAL is - begin - return Arccos (A, 360.0); - end ACOSD; - - function Arctan (A : REAL) return REAL is - begin - return Arctan (A, 1.0, 360.0); - end Arctan; - - function ATAND (A : REAL) return REAL is - begin - return Arctan (A, 1.0, 360.0); - end ATAND; - - function ATAN2D (A1, A2 : REAL) return REAL is - begin - return Arctan (A1, A2, 360.0); - end ATAN2D; - -end Math_Lib; |