------------------------------------------------------------------------------ -- -- -- GNAT RUNTIME COMPONENTS -- -- -- -- M A T H _ L I B -- -- -- -- B o d y -- -- -- -- Copyright (C) 1992-2000 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. -- 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 cyele 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;