diff options
Diffstat (limited to 'gcc-4.4.3/gcc/ada/a-ngelfu.adb')
| -rw-r--r-- | gcc-4.4.3/gcc/ada/a-ngelfu.adb | 1034 |
1 files changed, 0 insertions, 1034 deletions
diff --git a/gcc-4.4.3/gcc/ada/a-ngelfu.adb b/gcc-4.4.3/gcc/ada/a-ngelfu.adb deleted file mode 100644 index 62d88a48d..000000000 --- a/gcc-4.4.3/gcc/ada/a-ngelfu.adb +++ /dev/null @@ -1,1034 +0,0 @@ ------------------------------------------------------------------------------- --- -- --- GNAT RUN-TIME COMPONENTS -- --- -- --- ADA.NUMERICS.GENERIC_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. -- --- -- ------------------------------------------------------------------------------- - --- 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. - --- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, --- sinh, cosh, tanh from C library via math.h - -with Ada.Numerics.Aux; - -package body Ada.Numerics.Generic_Elementary_Functions is - - use type Ada.Numerics.Aux.Double; - - Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696; - Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755; - Half_Log_Two : constant := Log_Two / 2; - - subtype T is Float_Type'Base; - subtype Double is Aux.Double; - - Two_Pi : constant T := 2.0 * Pi; - Half_Pi : constant T := Pi / 2.0; - - Half_Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two; - Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two; - Sqrt_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa); - - ----------------------- - -- Local Subprograms -- - ----------------------- - - function Exp_Strict (X : Float_Type'Base) return Float_Type'Base; - -- Cody/Waite routine, supposedly more precise than the library - -- version. Currently only needed for Sinh/Cosh on X86 with the largest - -- FP type. - - function Local_Atan - (Y : Float_Type'Base; - X : Float_Type'Base := 1.0) - return Float_Type'Base; - -- Common code for arc tangent after cycle reduction - - ---------- - -- "**" -- - ---------- - - function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is - A_Right : Float_Type'Base; - Int_Part : Integer; - Result : Float_Type'Base; - R1 : Float_Type'Base; - Rest : Float_Type'Base; - - begin - if Left = 0.0 - and then Right = 0.0 - then - raise Argument_Error; - - elsif Left < 0.0 then - raise Argument_Error; - - elsif Right = 0.0 then - return 1.0; - - elsif Left = 0.0 then - if Right < 0.0 then - raise Constraint_Error; - else - return 0.0; - end if; - - elsif Left = 1.0 then - return 1.0; - - elsif Right = 1.0 then - return Left; - - else - begin - if Right = 2.0 then - return Left * Left; - - elsif Right = 0.5 then - return Sqrt (Left); - - else - A_Right := abs (Right); - - -- If exponent is larger than one, compute integer exponen- - -- tiation if possible, and evaluate fractional part with - -- more precision. The relative error is now proportional - -- to the fractional part of the exponent only. - - if A_Right > 1.0 - and then A_Right < Float_Type'Base (Integer'Last) - then - Int_Part := Integer (Float_Type'Base'Truncation (A_Right)); - Result := Left ** Int_Part; - Rest := A_Right - Float_Type'Base (Int_Part); - - -- Compute with two leading bits of the mantissa using - -- square roots. Bound to be better than logarithms, and - -- easily extended to greater precision. - - if Rest >= 0.5 then - R1 := Sqrt (Left); - Result := Result * R1; - Rest := Rest - 0.5; - - if Rest >= 0.25 then - Result := Result * Sqrt (R1); - Rest := Rest - 0.25; - end if; - - elsif Rest >= 0.25 then - Result := Result * Sqrt (Sqrt (Left)); - Rest := Rest - 0.25; - end if; - - Result := Result * - Float_Type'Base (Aux.Pow (Double (Left), Double (Rest))); - - if Right >= 0.0 then - return Result; - else - return (1.0 / Result); - end if; - else - return - Float_Type'Base (Aux.Pow (Double (Left), Double (Right))); - end if; - end if; - - exception - when others => - raise Constraint_Error; - end; - end if; - end "**"; - - ------------ - -- Arccos -- - ------------ - - -- Natural cycle - - function Arccos (X : Float_Type'Base) return Float_Type'Base is - Temp : Float_Type'Base; - - begin - if abs X > 1.0 then - raise Argument_Error; - - elsif abs X < Sqrt_Epsilon then - return Pi / 2.0 - X; - - elsif X = 1.0 then - return 0.0; - - elsif X = -1.0 then - return Pi; - end if; - - Temp := Float_Type'Base (Aux.Acos (Double (X))); - - if Temp < 0.0 then - Temp := Pi + Temp; - end if; - - return Temp; - end Arccos; - - -- Arbitrary cycle - - function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is - Temp : Float_Type'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif abs X > 1.0 then - raise Argument_Error; - - elsif abs X < Sqrt_Epsilon then - return Cycle / 4.0; - - elsif X = 1.0 then - return 0.0; - - elsif X = -1.0 then - return Cycle / 2.0; - end if; - - Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle); - - if Temp < 0.0 then - Temp := Cycle / 2.0 + Temp; - end if; - - return Temp; - end Arccos; - - ------------- - -- Arccosh -- - ------------- - - function Arccosh (X : Float_Type'Base) return Float_Type'Base is - begin - -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or - -- the proper approximation for X close to 1 or >> 1. - - if X < 1.0 then - raise Argument_Error; - - elsif X < 1.0 + Sqrt_Epsilon then - return Sqrt (2.0 * (X - 1.0)); - - elsif X > 1.0 / Sqrt_Epsilon then - return Log (X) + Log_Two; - - else - return Log (X + Sqrt ((X - 1.0) * (X + 1.0))); - end if; - end Arccosh; - - ------------ - -- Arccot -- - ------------ - - -- Natural cycle - - function Arccot - (X : Float_Type'Base; - Y : Float_Type'Base := 1.0) - return Float_Type'Base - is - begin - -- Just reverse arguments - - return Arctan (Y, X); - end Arccot; - - -- Arbitrary cycle - - function Arccot - (X : Float_Type'Base; - Y : Float_Type'Base := 1.0; - Cycle : Float_Type'Base) - return Float_Type'Base - is - begin - -- Just reverse arguments - - return Arctan (Y, X, Cycle); - end Arccot; - - ------------- - -- Arccoth -- - ------------- - - function Arccoth (X : Float_Type'Base) return Float_Type'Base is - begin - if abs X > 2.0 then - return Arctanh (1.0 / X); - - elsif abs X = 1.0 then - raise Constraint_Error; - - elsif abs X < 1.0 then - raise Argument_Error; - - else - -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the - -- other has error 0 or Epsilon. - - return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0))); - end if; - end Arccoth; - - ------------ - -- Arcsin -- - ------------ - - -- Natural cycle - - function Arcsin (X : Float_Type'Base) return Float_Type'Base is - begin - if abs X > 1.0 then - raise Argument_Error; - - elsif abs X < Sqrt_Epsilon then - return X; - - elsif X = 1.0 then - return Pi / 2.0; - - elsif X = -1.0 then - return -(Pi / 2.0); - end if; - - return Float_Type'Base (Aux.Asin (Double (X))); - end Arcsin; - - -- Arbitrary cycle - - function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif abs X > 1.0 then - raise Argument_Error; - - elsif X = 0.0 then - return X; - - elsif X = 1.0 then - return Cycle / 4.0; - - elsif X = -1.0 then - return -(Cycle / 4.0); - end if; - - return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle); - end Arcsin; - - ------------- - -- Arcsinh -- - ------------- - - function Arcsinh (X : Float_Type'Base) return Float_Type'Base is - begin - if abs X < Sqrt_Epsilon then - return X; - - elsif X > 1.0 / Sqrt_Epsilon then - return Log (X) + Log_Two; - - elsif X < -(1.0 / Sqrt_Epsilon) then - return -(Log (-X) + Log_Two); - - elsif X < 0.0 then - return -Log (abs X + Sqrt (X * X + 1.0)); - - else - return Log (X + Sqrt (X * X + 1.0)); - end if; - end Arcsinh; - - ------------ - -- Arctan -- - ------------ - - -- Natural cycle - - function Arctan - (Y : Float_Type'Base; - X : Float_Type'Base := 1.0) - return Float_Type'Base - is - begin - if X = 0.0 - and then Y = 0.0 - then - raise Argument_Error; - - elsif Y = 0.0 then - if X > 0.0 then - return 0.0; - else -- X < 0.0 - return Pi * Float_Type'Copy_Sign (1.0, Y); - end if; - - elsif X = 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, X); - end if; - end Arctan; - - -- Arbitrary cycle - - function Arctan - (Y : Float_Type'Base; - X : Float_Type'Base := 1.0; - Cycle : Float_Type'Base) - return Float_Type'Base - is - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif X = 0.0 - and then Y = 0.0 - then - raise Argument_Error; - - elsif Y = 0.0 then - if X > 0.0 then - return 0.0; - else -- X < 0.0 - return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y); - end if; - - elsif X = 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, X) * Cycle / Two_Pi; - end if; - end Arctan; - - ------------- - -- Arctanh -- - ------------- - - function Arctanh (X : Float_Type'Base) return Float_Type'Base is - A, B, D, A_Plus_1, A_From_1 : Float_Type'Base; - Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa; - - begin - -- The naive formula: - - -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X) - - -- is not well-behaved numerically when X < 0.5 and when X is close - -- to one. The following is accurate but probably not optimal. - - if abs X = 1.0 then - raise Constraint_Error; - - elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then - - if abs X >= 1.0 then - raise Argument_Error; - else - - -- The one case that overflows if put through the method below: - -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is - -- accurate. This simplifies to: - - return Float_Type'Copy_Sign ( - Half_Log_Two * Float_Type'Base (Mantissa + 1), X); - end if; - - -- elsif abs X <= 0.5 then - -- why is above line commented out ??? - - else - -- Use several piecewise linear approximations. - -- A is close to X, chosen so 1.0 + A, 1.0 - A, and X - A are exact. - -- The two scalings remove the low-order bits of X. - - A := Float_Type'Base'Scaling ( - Float_Type'Base (Long_Long_Integer - (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa); - - B := X - A; -- This is exact; abs B <= 2**(-Mantissa). - A_Plus_1 := 1.0 + A; -- This is exact. - A_From_1 := 1.0 - A; -- Ditto. - D := A_Plus_1 * A_From_1; -- 1 - A*A. - - -- use one term of the series expansion: - -- f (x + e) = f(x) + e * f'(x) + .. - - -- The derivative of Arctanh at A is 1/(1-A*A). Next term is - -- A*(B/D)**2 (if a quadratic approximation is ever needed). - - return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D; - - -- else - -- return 0.5 * Log ((X + 1.0) / (1.0 - X)); - -- why are above lines commented out ??? - end if; - end Arctanh; - - --------- - -- Cos -- - --------- - - -- Natural cycle - - function Cos (X : Float_Type'Base) return Float_Type'Base is - begin - if X = 0.0 then - return 1.0; - - elsif abs X < Sqrt_Epsilon then - return 1.0; - - end if; - - return Float_Type'Base (Aux.Cos (Double (X))); - end Cos; - - -- Arbitrary cycle - - function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is - begin - -- Just reuse the code for Sin. The potential small - -- loss of speed is negligible with proper (front-end) inlining. - - return -Sin (abs X - Cycle * 0.25, Cycle); - end Cos; - - ---------- - -- Cosh -- - ---------- - - function Cosh (X : Float_Type'Base) return Float_Type'Base is - Lnv : constant Float_Type'Base := 8#0.542714#; - V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4; - Y : constant Float_Type'Base := abs X; - Z : Float_Type'Base; - - begin - if Y < Sqrt_Epsilon then - return 1.0; - - elsif Y > Log_Inverse_Epsilon then - Z := Exp_Strict (Y - Lnv); - return (Z + V2minus1 * Z); - - else - Z := Exp_Strict (Y); - return 0.5 * (Z + 1.0 / Z); - end if; - - end Cosh; - - --------- - -- Cot -- - --------- - - -- Natural cycle - - function Cot (X : Float_Type'Base) return Float_Type'Base is - begin - if X = 0.0 then - raise Constraint_Error; - - elsif abs X < Sqrt_Epsilon then - return 1.0 / X; - end if; - - return 1.0 / Float_Type'Base (Aux.Tan (Double (X))); - end Cot; - - -- Arbitrary cycle - - function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is - T : Float_Type'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - end if; - - T := Float_Type'Base'Remainder (X, Cycle); - - if T = 0.0 or abs T = 0.5 * Cycle then - raise Constraint_Error; - - elsif abs T < Sqrt_Epsilon then - return 1.0 / T; - - elsif abs T = 0.25 * Cycle then - return 0.0; - - else - T := T / Cycle * Two_Pi; - return Cos (T) / Sin (T); - end if; - end Cot; - - ---------- - -- Coth -- - ---------- - - function Coth (X : Float_Type'Base) return Float_Type'Base is - begin - if X = 0.0 then - raise Constraint_Error; - - elsif X < Half_Log_Epsilon then - return -1.0; - - elsif X > -Half_Log_Epsilon then - return 1.0; - - elsif abs X < Sqrt_Epsilon then - return 1.0 / X; - end if; - - return 1.0 / Float_Type'Base (Aux.Tanh (Double (X))); - end Coth; - - --------- - -- Exp -- - --------- - - function Exp (X : Float_Type'Base) return Float_Type'Base is - Result : Float_Type'Base; - - begin - if X = 0.0 then - return 1.0; - end if; - - Result := Float_Type'Base (Aux.Exp (Double (X))); - - -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows - -- is False, then we can just leave it as an infinity (and indeed we - -- prefer to do so). But if Machine_Overflows is True, then we have - -- to raise a Constraint_Error exception as required by the RM. - - if Float_Type'Machine_Overflows and then not Result'Valid then - raise Constraint_Error; - end if; - - return Result; - end Exp; - - ---------------- - -- Exp_Strict -- - ---------------- - - function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is - G : Float_Type'Base; - Z : Float_Type'Base; - - P0 : constant := 0.25000_00000_00000_00000; - P1 : constant := 0.75753_18015_94227_76666E-2; - P2 : constant := 0.31555_19276_56846_46356E-4; - - Q0 : constant := 0.5; - Q1 : constant := 0.56817_30269_85512_21787E-1; - Q2 : constant := 0.63121_89437_43985_02557E-3; - Q3 : constant := 0.75104_02839_98700_46114E-6; - - C1 : constant := 8#0.543#; - C2 : constant := -2.1219_44400_54690_58277E-4; - Le : constant := 1.4426_95040_88896_34074; - - XN : Float_Type'Base; - P, Q, R : Float_Type'Base; - - begin - if X = 0.0 then - return 1.0; - end if; - - XN := Float_Type'Base'Rounding (X * Le); - G := (X - XN * C1) - XN * C2; - Z := G * G; - P := G * ((P2 * Z + P1) * Z + P0); - Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0; - R := 0.5 + P / (Q - P); - - R := Float_Type'Base'Scaling (R, Integer (XN) + 1); - - -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows - -- is False, then we can just leave it as an infinity (and indeed we - -- prefer to do so). But if Machine_Overflows is True, then we have - -- to raise a Constraint_Error exception as required by the RM. - - if Float_Type'Machine_Overflows and then not R'Valid then - raise Constraint_Error; - else - return R; - end if; - - end Exp_Strict; - - ---------------- - -- Local_Atan -- - ---------------- - - function Local_Atan - (Y : Float_Type'Base; - X : Float_Type'Base := 1.0) - return Float_Type'Base - is - Z : Float_Type'Base; - Raw_Atan : Float_Type'Base; - - begin - if abs Y > abs X then - Z := abs (X / Y); - else - Z := abs (Y / X); - end if; - - if Z < Sqrt_Epsilon then - Raw_Atan := Z; - - elsif Z = 1.0 then - Raw_Atan := Pi / 4.0; - - else - Raw_Atan := Float_Type'Base (Aux.Atan (Double (Z))); - end if; - - if abs Y > abs X then - Raw_Atan := Half_Pi - Raw_Atan; - end if; - - if X > 0.0 then - if Y > 0.0 then - return Raw_Atan; - else -- Y < 0.0 - return -Raw_Atan; - end if; - - else -- X < 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 (X : Float_Type'Base) return Float_Type'Base is - begin - if X < 0.0 then - raise Argument_Error; - - elsif X = 0.0 then - raise Constraint_Error; - - elsif X = 1.0 then - return 0.0; - end if; - - return Float_Type'Base (Aux.Log (Double (X))); - end Log; - - -- Arbitrary base - - function Log (X, Base : Float_Type'Base) return Float_Type'Base is - begin - if X < 0.0 then - raise Argument_Error; - - elsif Base <= 0.0 or else Base = 1.0 then - raise Argument_Error; - - elsif X = 0.0 then - raise Constraint_Error; - - elsif X = 1.0 then - return 0.0; - end if; - - return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base))); - end Log; - - --------- - -- Sin -- - --------- - - -- Natural cycle - - function Sin (X : Float_Type'Base) return Float_Type'Base is - begin - if abs X < Sqrt_Epsilon then - return X; - end if; - - return Float_Type'Base (Aux.Sin (Double (X))); - end Sin; - - -- Arbitrary cycle - - function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is - T : Float_Type'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif X = 0.0 then - -- Is this test really needed on any machine ??? - return X; - end if; - - T := Float_Type'Base'Remainder (X, Cycle); - - -- The following two reductions reduce the argument - -- to the interval [-0.25 * Cycle, 0.25 * Cycle]. - -- This reduction is exact and is needed to prevent - -- inaccuracy that may result if the sinus function - -- a different (more accurate) value of Pi in its - -- reduction than is used in the multiplication with Two_Pi. - - if abs T > 0.25 * Cycle then - T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T; - end if; - - -- Could test for 12.0 * abs T = Cycle, and return - -- an exact value in those cases. It is not clear that - -- this is worth the extra test though. - - return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi))); - end Sin; - - ---------- - -- Sinh -- - ---------- - - function Sinh (X : Float_Type'Base) return Float_Type'Base is - Lnv : constant Float_Type'Base := 8#0.542714#; - V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4; - Y : constant Float_Type'Base := abs X; - F : constant Float_Type'Base := Y * Y; - Z : Float_Type'Base; - - Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7; - - begin - if Y < Sqrt_Epsilon then - return X; - - elsif Y > Log_Inverse_Epsilon then - Z := Exp_Strict (Y - Lnv); - Z := Z + V2minus1 * Z; - - elsif Y < 1.0 then - - if Float_Digits_1_6 then - - -- Use expansion provided by Cody and Waite, p. 226. Note that - -- leading term of the polynomial in Q is exactly 1.0. - - declare - P0 : constant := -0.71379_3159E+1; - P1 : constant := -0.19033_3399E+0; - Q0 : constant := -0.42827_7109E+2; - - begin - Z := Y + Y * F * (P1 * F + P0) / (F + Q0); - end; - - else - declare - P0 : constant := -0.35181_28343_01771_17881E+6; - P1 : constant := -0.11563_52119_68517_68270E+5; - P2 : constant := -0.16375_79820_26307_51372E+3; - P3 : constant := -0.78966_12741_73570_99479E+0; - Q0 : constant := -0.21108_77005_81062_71242E+7; - Q1 : constant := 0.36162_72310_94218_36460E+5; - Q2 : constant := -0.27773_52311_96507_01667E+3; - - begin - Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0) - / (((F + Q2) * F + Q1) * F + Q0); - end; - end if; - - else - Z := Exp_Strict (Y); - Z := 0.5 * (Z - 1.0 / Z); - end if; - - if X > 0.0 then - return Z; - else - return -Z; - end if; - end Sinh; - - ---------- - -- Sqrt -- - ---------- - - function Sqrt (X : Float_Type'Base) return Float_Type'Base is - begin - if X < 0.0 then - raise Argument_Error; - - -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE - - elsif X = 0.0 then - return X; - - end if; - - return Float_Type'Base (Aux.Sqrt (Double (X))); - end Sqrt; - - --------- - -- Tan -- - --------- - - -- Natural cycle - - function Tan (X : Float_Type'Base) return Float_Type'Base is - begin - if abs X < Sqrt_Epsilon then - return X; - - elsif abs X = Pi / 2.0 then - raise Constraint_Error; - end if; - - return Float_Type'Base (Aux.Tan (Double (X))); - end Tan; - - -- Arbitrary cycle - - function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is - T : Float_Type'Base; - - begin - if Cycle <= 0.0 then - raise Argument_Error; - - elsif X = 0.0 then - return X; - end if; - - T := Float_Type'Base'Remainder (X, Cycle); - - if abs T = 0.25 * Cycle then - raise Constraint_Error; - - elsif abs T = 0.5 * Cycle then - return 0.0; - - else - T := T / Cycle * Two_Pi; - return Sin (T) / Cos (T); - end if; - - end Tan; - - ---------- - -- Tanh -- - ---------- - - function Tanh (X : Float_Type'Base) return Float_Type'Base is - P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4; - P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2; - P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0; - - Q0 : constant Float_Type'Base := 0.48402_35707_19886_88686E+4; - Q1 : constant Float_Type'Base := 0.22337_72071_89623_12926E+4; - Q2 : constant Float_Type'Base := 0.11274_47438_05349_49335E+3; - Q3 : constant Float_Type'Base := 0.10000_00000_00000_00000E+1; - - Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570; - - P, Q, R : Float_Type'Base; - Y : constant Float_Type'Base := abs X; - G : constant Float_Type'Base := Y * Y; - - Float_Type_Digits_15_Or_More : constant Boolean := - Float_Type'Digits > 14; - - begin - if X < Half_Log_Epsilon then - return -1.0; - - elsif X > -Half_Log_Epsilon then - return 1.0; - - elsif Y < Sqrt_Epsilon then - return X; - - elsif Y < Half_Ln3 - and then Float_Type_Digits_15_Or_More - then - P := (P2 * G + P1) * G + P0; - Q := ((Q3 * G + Q2) * G + Q1) * G + Q0; - R := G * (P / Q); - return X + X * R; - - else - return Float_Type'Base (Aux.Tanh (Double (X))); - end if; - end Tanh; - -end Ada.Numerics.Generic_Elementary_Functions; |