------------------------------------------------------------------------------ -- -- -- GNAT RUNTIME COMPONENTS -- -- -- -- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS -- -- -- -- B o d y -- -- -- -- Copyright (C) 1992-2004, 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 2, 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. See the GNU General Public License -- -- for more details. You should have received a copy of the GNU General -- -- Public License distributed with GNAT; see file COPYING. If not, write -- -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, -- -- MA 02111-1307, USA. -- -- -- -- As a special exception, if other files instantiate generics from this -- -- unit, or you link this unit with other files to produce an executable, -- -- this unit does not by itself cause the resulting executable to be -- -- covered by the GNU General Public License. This exception does not -- -- however invalidate any other reasons why the executable file might be -- -- covered by the GNU Public License. -- -- -- -- 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 cyele 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;