565 lines
18 KiB
Ada
565 lines
18 KiB
Ada
|
------------------------------------------------------------------------------
|
||
|
|
-- --
|
||
|
|
-- GNAT RUN-TIME COMPONENTS --
|
||
|
|
-- --
|
||
|
|
-- A D A . N U M E R I C S . A U X --
|
||
|
|
-- --
|
||
|
|
-- B o d y --
|
||
|
|
-- (Machine Version for x86) --
|
||
|
|
-- --
|
||
|
|
-- Copyright (C) 1998-2014, 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. --
|
||
|
|
-- --
|
||
|
|
------------------------------------------------------------------------------
|
||
|
|
|
||
|
|
with System.Machine_Code; use System.Machine_Code;
|
||
|
|
|
||
|
|
package body Ada.Numerics.Aux is
|
||
|
|
|
||
|
|
NL : constant String := ASCII.LF & ASCII.HT;
|
||
|
|
|
||
|
|
-----------------------
|
||
|
|
-- Local subprograms --
|
||
|
|
-----------------------
|
||
|
|
|
||
|
|
function Is_Nan (X : Double) return Boolean;
|
||
|
|
-- Return True iff X is a IEEE NaN value
|
||
|
|
|
||
|
|
function Logarithmic_Pow (X, Y : Double) return Double;
|
||
|
|
-- Implementation of X**Y using Exp and Log functions (binary base)
|
||
|
|
-- to calculate the exponentiation. This is used by Pow for values
|
||
|
|
-- for values of Y in the open interval (-0.25, 0.25)
|
||
|
|
|
||
|
|
procedure Reduce (X : in out Double; Q : out Natural);
|
||
|
|
-- Implements reduction of X by Pi/2. Q is the quadrant of the final
|
||
|
|
-- result in the range 0 .. 3. The absolute value of X is at most Pi.
|
||
|
|
|
||
|
|
pragma Inline (Is_Nan);
|
||
|
|
pragma Inline (Reduce);
|
||
|
|
|
||
|
|
--------------------------------
|
||
|
|
-- Basic Elementary Functions --
|
||
|
|
--------------------------------
|
||
|
|
|
||
|
|
-- This section implements a few elementary functions that are used to
|
||
|
|
-- build the more complex ones. This ordering enables better inlining.
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Atan --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Atan (X : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
|
||
|
|
begin
|
||
|
|
Asm (Template =>
|
||
|
|
"fld1" & NL
|
||
|
|
& "fpatan",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", X));
|
||
|
|
|
||
|
|
-- The result value is NaN iff input was invalid
|
||
|
|
|
||
|
|
if not (Result = Result) then
|
||
|
|
raise Argument_Error;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Atan;
|
||
|
|
|
||
|
|
---------
|
||
|
|
-- Exp --
|
||
|
|
---------
|
||
|
|
|
||
|
|
function Exp (X : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
begin
|
||
|
|
Asm (Template =>
|
||
|
|
"fldl2e " & NL
|
||
|
|
& "fmulp %%st, %%st(1)" & NL -- X * log2 (E)
|
||
|
|
& "fld %%st(0) " & NL
|
||
|
|
& "frndint " & NL -- Integer (X * Log2 (E))
|
||
|
|
& "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
|
||
|
|
& "fxch " & NL
|
||
|
|
& "f2xm1 " & NL -- 2**(...) - 1
|
||
|
|
& "fld1 " & NL
|
||
|
|
& "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
|
||
|
|
& "fscale " & NL -- E ** X
|
||
|
|
& "fstp %%st(1) ",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", X));
|
||
|
|
return Result;
|
||
|
|
end Exp;
|
||
|
|
|
||
|
|
------------
|
||
|
|
-- Is_Nan --
|
||
|
|
------------
|
||
|
|
|
||
|
|
function Is_Nan (X : Double) return Boolean is
|
||
|
|
begin
|
||
|
|
-- The IEEE NaN values are the only ones that do not equal themselves
|
||
|
|
|
||
|
|
return not (X = X);
|
||
|
|
end Is_Nan;
|
||
|
|
|
||
|
|
---------
|
||
|
|
-- Log --
|
||
|
|
---------
|
||
|
|
|
||
|
|
function Log (X : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
|
||
|
|
begin
|
||
|
|
Asm (Template =>
|
||
|
|
"fldln2 " & NL
|
||
|
|
& "fxch " & NL
|
||
|
|
& "fyl2x " & NL,
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", X));
|
||
|
|
return Result;
|
||
|
|
end Log;
|
||
|
|
|
||
|
|
------------
|
||
|
|
-- Reduce --
|
||
|
|
------------
|
||
|
|
|
||
|
|
procedure Reduce (X : in out Double; Q : out Natural) is
|
||
|
|
Half_Pi : constant := Pi / 2.0;
|
||
|
|
Two_Over_Pi : constant := 2.0 / Pi;
|
||
|
|
|
||
|
|
HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
|
||
|
|
M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
|
||
|
|
P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
|
||
|
|
P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
|
||
|
|
P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
|
||
|
|
P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
|
||
|
|
P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
|
||
|
|
- P4, HM);
|
||
|
|
P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
|
||
|
|
K : Double := X * Two_Over_Pi;
|
||
|
|
begin
|
||
|
|
-- For X < 2.0**32, all products below are computed exactly.
|
||
|
|
-- Due to cancellation effects all subtractions are exact as well.
|
||
|
|
-- As no double extended floating-point number has more than 75
|
||
|
|
-- zeros after the binary point, the result will be the correctly
|
||
|
|
-- rounded result of X - K * (Pi / 2.0).
|
||
|
|
|
||
|
|
while abs K >= 2.0**HM loop
|
||
|
|
K := K * M - (K * M - K);
|
||
|
|
X := (((((X - K * P1) - K * P2) - K * P3)
|
||
|
|
- K * P4) - K * P5) - K * P6;
|
||
|
|
K := X * Two_Over_Pi;
|
||
|
|
end loop;
|
||
|
|
|
||
|
|
if K /= K then
|
||
|
|
|
||
|
|
-- K is not a number, because X was not finite
|
||
|
|
|
||
|
|
raise Constraint_Error;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
K := Double'Rounding (K);
|
||
|
|
Q := Integer (K) mod 4;
|
||
|
|
X := (((((X - K * P1) - K * P2) - K * P3)
|
||
|
|
- K * P4) - K * P5) - K * P6;
|
||
|
|
end Reduce;
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Sqrt --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Sqrt (X : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
|
||
|
|
begin
|
||
|
|
if X < 0.0 then
|
||
|
|
raise Argument_Error;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
Asm (Template => "fsqrt",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", X));
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Sqrt;
|
||
|
|
|
||
|
|
--------------------------------
|
||
|
|
-- Other Elementary Functions --
|
||
|
|
--------------------------------
|
||
|
|
|
||
|
|
-- These are built using the previously implemented basic functions
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Acos --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Acos (X : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
|
||
|
|
begin
|
||
|
|
Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
|
||
|
|
|
||
|
|
-- The result value is NaN iff input was invalid
|
||
|
|
|
||
|
|
if Is_Nan (Result) then
|
||
|
|
raise Argument_Error;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Acos;
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Asin --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Asin (X : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
|
||
|
|
begin
|
||
|
|
Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
|
||
|
|
|
||
|
|
-- The result value is NaN iff input was invalid
|
||
|
|
|
||
|
|
if Is_Nan (Result) then
|
||
|
|
raise Argument_Error;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Asin;
|
||
|
|
|
||
|
|
---------
|
||
|
|
-- Cos --
|
||
|
|
---------
|
||
|
|
|
||
|
|
function Cos (X : Double) return Double is
|
||
|
|
Reduced_X : Double := abs X;
|
||
|
|
Result : Double;
|
||
|
|
Quadrant : Natural range 0 .. 3;
|
||
|
|
|
||
|
|
begin
|
||
|
|
if Reduced_X > Pi / 4.0 then
|
||
|
|
Reduce (Reduced_X, Quadrant);
|
||
|
|
|
||
|
|
case Quadrant is
|
||
|
|
when 0 =>
|
||
|
|
Asm (Template => "fcos",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
when 1 =>
|
||
|
|
Asm (Template => "fsin",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", -Reduced_X));
|
||
|
|
when 2 =>
|
||
|
|
Asm (Template => "fcos ; fchs",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
when 3 =>
|
||
|
|
Asm (Template => "fsin",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
end case;
|
||
|
|
|
||
|
|
else
|
||
|
|
Asm (Template => "fcos",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Cos;
|
||
|
|
|
||
|
|
---------------------
|
||
|
|
-- Logarithmic_Pow --
|
||
|
|
---------------------
|
||
|
|
|
||
|
|
function Logarithmic_Pow (X, Y : Double) return Double is
|
||
|
|
Result : Double;
|
||
|
|
begin
|
||
|
|
Asm (Template => "" -- X : Y
|
||
|
|
& "fyl2x " & NL -- Y * Log2 (X)
|
||
|
|
& "fld %%st(0) " & NL -- Y * Log2 (X) : Y * Log2 (X)
|
||
|
|
& "frndint " & NL -- Int (...) : Y * Log2 (X)
|
||
|
|
& "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...)
|
||
|
|
& "fxch " & NL -- Fract (...) : Int (...)
|
||
|
|
& "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...)
|
||
|
|
& "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...)
|
||
|
|
& "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...)
|
||
|
|
& "fscale ", -- 2**(Fract (...) + Int (...))
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs =>
|
||
|
|
(Double'Asm_Input ("0", X),
|
||
|
|
Double'Asm_Input ("u", Y)));
|
||
|
|
return Result;
|
||
|
|
end Logarithmic_Pow;
|
||
|
|
|
||
|
|
---------
|
||
|
|
-- Pow --
|
||
|
|
---------
|
||
|
|
|
||
|
|
function Pow (X, Y : Double) return Double is
|
||
|
|
type Mantissa_Type is mod 2**Double'Machine_Mantissa;
|
||
|
|
-- Modular type that can hold all bits of the mantissa of Double
|
||
|
|
|
||
|
|
-- For negative exponents, do divide at the end of the processing
|
||
|
|
|
||
|
|
Negative_Y : constant Boolean := Y < 0.0;
|
||
|
|
Abs_Y : constant Double := abs Y;
|
||
|
|
|
||
|
|
-- During this function the following invariant is kept:
|
||
|
|
-- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
|
||
|
|
|
||
|
|
Base : Double := X;
|
||
|
|
|
||
|
|
Exp_High : Double := Double'Floor (Abs_Y);
|
||
|
|
Exp_Mid : Double;
|
||
|
|
Exp_Low : Double;
|
||
|
|
Exp_Int : Mantissa_Type;
|
||
|
|
|
||
|
|
Factor : Double := 1.0;
|
||
|
|
|
||
|
|
begin
|
||
|
|
-- Select algorithm for calculating Pow (integer cases fall through)
|
||
|
|
|
||
|
|
if Exp_High >= 2.0**Double'Machine_Mantissa then
|
||
|
|
|
||
|
|
-- In case of Y that is IEEE infinity, just raise constraint error
|
||
|
|
|
||
|
|
if Exp_High > Double'Safe_Last then
|
||
|
|
raise Constraint_Error;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
-- Large values of Y are even integers and will stay integer
|
||
|
|
-- after division by two.
|
||
|
|
|
||
|
|
loop
|
||
|
|
-- Exp_Mid and Exp_Low are zero, so
|
||
|
|
-- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
|
||
|
|
|
||
|
|
Exp_High := Exp_High / 2.0;
|
||
|
|
Base := Base * Base;
|
||
|
|
exit when Exp_High < 2.0**Double'Machine_Mantissa;
|
||
|
|
end loop;
|
||
|
|
|
||
|
|
elsif Exp_High /= Abs_Y then
|
||
|
|
Exp_Low := Abs_Y - Exp_High;
|
||
|
|
Factor := 1.0;
|
||
|
|
|
||
|
|
if Exp_Low /= 0.0 then
|
||
|
|
|
||
|
|
-- Exp_Low now is in interval (0.0, 1.0)
|
||
|
|
-- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
|
||
|
|
|
||
|
|
Exp_Mid := 0.0;
|
||
|
|
Exp_Low := Exp_Low - Exp_Mid;
|
||
|
|
|
||
|
|
if Exp_Low >= 0.5 then
|
||
|
|
Factor := Sqrt (X);
|
||
|
|
Exp_Low := Exp_Low - 0.5; -- exact
|
||
|
|
|
||
|
|
if Exp_Low >= 0.25 then
|
||
|
|
Factor := Factor * Sqrt (Factor);
|
||
|
|
Exp_Low := Exp_Low - 0.25; -- exact
|
||
|
|
end if;
|
||
|
|
|
||
|
|
elsif Exp_Low >= 0.25 then
|
||
|
|
Factor := Sqrt (Sqrt (X));
|
||
|
|
Exp_Low := Exp_Low - 0.25; -- exact
|
||
|
|
end if;
|
||
|
|
|
||
|
|
-- Exp_Low now is in interval (0.0, 0.25)
|
||
|
|
|
||
|
|
-- This means it is safe to call Logarithmic_Pow
|
||
|
|
-- for the remaining part.
|
||
|
|
|
||
|
|
Factor := Factor * Logarithmic_Pow (X, Exp_Low);
|
||
|
|
end if;
|
||
|
|
|
||
|
|
elsif X = 0.0 then
|
||
|
|
return 0.0;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
-- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
|
||
|
|
|
||
|
|
Exp_Int := Mantissa_Type (Exp_High);
|
||
|
|
|
||
|
|
-- Standard way for processing integer powers > 0
|
||
|
|
|
||
|
|
while Exp_Int > 1 loop
|
||
|
|
if (Exp_Int and 1) = 1 then
|
||
|
|
|
||
|
|
-- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
|
||
|
|
|
||
|
|
Factor := Factor * Base;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
-- Exp_Int is even and Exp_Int > 0, so
|
||
|
|
-- Base**Y = (Base**2)**(Exp_Int / 2)
|
||
|
|
|
||
|
|
Base := Base * Base;
|
||
|
|
Exp_Int := Exp_Int / 2;
|
||
|
|
end loop;
|
||
|
|
|
||
|
|
-- Exp_Int = 1 or Exp_Int = 0
|
||
|
|
|
||
|
|
if Exp_Int = 1 then
|
||
|
|
Factor := Base * Factor;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
if Negative_Y then
|
||
|
|
Factor := 1.0 / Factor;
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Factor;
|
||
|
|
end Pow;
|
||
|
|
|
||
|
|
---------
|
||
|
|
-- Sin --
|
||
|
|
---------
|
||
|
|
|
||
|
|
function Sin (X : Double) return Double is
|
||
|
|
Reduced_X : Double := X;
|
||
|
|
Result : Double;
|
||
|
|
Quadrant : Natural range 0 .. 3;
|
||
|
|
|
||
|
|
begin
|
||
|
|
if abs X > Pi / 4.0 then
|
||
|
|
Reduce (Reduced_X, Quadrant);
|
||
|
|
|
||
|
|
case Quadrant is
|
||
|
|
when 0 =>
|
||
|
|
Asm (Template => "fsin",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
when 1 =>
|
||
|
|
Asm (Template => "fcos",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
when 2 =>
|
||
|
|
Asm (Template => "fsin",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", -Reduced_X));
|
||
|
|
when 3 =>
|
||
|
|
Asm (Template => "fcos ; fchs",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
end case;
|
||
|
|
|
||
|
|
else
|
||
|
|
Asm (Template => "fsin",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Sin;
|
||
|
|
|
||
|
|
---------
|
||
|
|
-- Tan --
|
||
|
|
---------
|
||
|
|
|
||
|
|
function Tan (X : Double) return Double is
|
||
|
|
Reduced_X : Double := X;
|
||
|
|
Result : Double;
|
||
|
|
Quadrant : Natural range 0 .. 3;
|
||
|
|
|
||
|
|
begin
|
||
|
|
if abs X > Pi / 4.0 then
|
||
|
|
Reduce (Reduced_X, Quadrant);
|
||
|
|
|
||
|
|
if Quadrant mod 2 = 0 then
|
||
|
|
Asm (Template => "fptan" & NL
|
||
|
|
& "ffree %%st(0)" & NL
|
||
|
|
& "fincstp",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
else
|
||
|
|
Asm (Template => "fsincos" & NL
|
||
|
|
& "fdivp %%st, %%st(1)" & NL
|
||
|
|
& "fchs",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
end if;
|
||
|
|
|
||
|
|
else
|
||
|
|
Asm (Template =>
|
||
|
|
"fptan " & NL
|
||
|
|
& "ffree %%st(0) " & NL
|
||
|
|
& "fincstp ",
|
||
|
|
Outputs => Double'Asm_Output ("=t", Result),
|
||
|
|
Inputs => Double'Asm_Input ("0", Reduced_X));
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return Result;
|
||
|
|
end Tan;
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Sinh --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Sinh (X : Double) return Double is
|
||
|
|
begin
|
||
|
|
-- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
|
||
|
|
|
||
|
|
if abs X < 25.0 then
|
||
|
|
return (Exp (X) - Exp (-X)) / 2.0;
|
||
|
|
else
|
||
|
|
return Exp (X) / 2.0;
|
||
|
|
end if;
|
||
|
|
end Sinh;
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Cosh --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Cosh (X : Double) return Double is
|
||
|
|
begin
|
||
|
|
-- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
|
||
|
|
|
||
|
|
if abs X < 22.0 then
|
||
|
|
return (Exp (X) + Exp (-X)) / 2.0;
|
||
|
|
else
|
||
|
|
return Exp (X) / 2.0;
|
||
|
|
end if;
|
||
|
|
end Cosh;
|
||
|
|
|
||
|
|
----------
|
||
|
|
-- Tanh --
|
||
|
|
----------
|
||
|
|
|
||
|
|
function Tanh (X : Double) return Double is
|
||
|
|
begin
|
||
|
|
-- Return the Hyperbolic Tangent of x
|
||
|
|
|
||
|
|
-- x -x
|
||
|
|
-- e - e Sinh (X)
|
||
|
|
-- Tanh (X) is defined to be ----------- = --------
|
||
|
|
-- x -x Cosh (X)
|
||
|
|
-- e + e
|
||
|
|
|
||
|
|
if abs X > 23.0 then
|
||
|
|
return Double'Copy_Sign (1.0, X);
|
||
|
|
end if;
|
||
|
|
|
||
|
|
return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
|
||
|
|
end Tanh;
|
||
|
|
|
||
|
|
end Ada.Numerics.Aux;
|