1028 lines
22 KiB
Ada
1028 lines
22 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT RUN-TIME COMPONENTS --
|
|
-- --
|
|
-- M A T H _ L I B --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 1992-2005 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, 51 Franklin Street, Fifth Floor, --
|
|
-- Boston, MA 02110-1301, 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;
|