+++ /dev/null
-------------------------------------------------------------------------------
--- --
--- 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;
projects / gcc.git / commitdiff