separate ( Generic_Elementary_Functions )

function Tanh( X : Float_Type ) return Float_Type is

-- On input, X is a floating-point value in Float_Type;
-- On output, the value of tanh(X) (the hyperbolic tangent of X) is returned.

-- The definition of tanh(Y) is sinh(Y)/cosh(Y), which is also equivalent
-- to the following three formulas.
--      1.  ( exp(Y) - exp(-Y) ) / ( exp(Y) + exp(-Y) )
--      2.  ( 1 - ( 2 / ( exp(2*Y) + 1 ) ) )
--      3.  ( exp(2*Y) - 1 ) / ( exp(2*Y) + 1 ).
-- but computationally, some formulas are better on some ranges.

   Z, Sign_Y                            : Common_Float;

   Y, Abs_Y                             : Common_Float;

   Log2         : constant Common_Float :=
                  16#0.B17217F7D1CF79ABC9E3B39803F2F6AF40#;

   Base_Digits  : constant Common_Float :=
                  Common_Float( 6 * Float_Type'Base'Digits );

   Log2_Times_2 : constant Common_Float := ( 2.0 * Log2 );

   Cond         : constant Common_Float := ( Base_Digits * Log2 );


begin

-- Filter out exceptional cases.

   if (X = 0.0) then
      return( X );
   end if;

   Y     := Common_Float( X );
   Abs_Y := abs(Y);

   if ( Y >= 0.0 ) then
      Sign_Y := 1.0;
   else
      Sign_Y := -1.0;
   end if;


   if ( Abs_Y <= ( Log2_Times_2 ) ) then
--    Formula 3 should be used in this situation to guarantee accuracy.

      Z      := KF_Em1 ( 2.0 * Abs_Y );
      Z      := Sign_Y * ( Z / ( Z + 2.0 ) );
      return ( Float_Type(Z) );

   elsif (Abs_Y > Cond) then
--    Formula 2 should be used in this situation to guarantee accuracy,
--    but observe that 2/(exp(2*Y) + 1) will be so small compared to 1
--    that it is negligible.

      return ( Float_Type(Sign_Y) );

   else
--   When ( Log2_Times_2 < Abs_Y <= Cond ), use formula 2 for best accuracy. 

      Z := KF_Em1 ( 2.0 * Abs_Y );
      Z := Sign_Y * ( 1.0 - 2.0 / ( Z + 2.0 ) );
      return ( Float_Type(Z) );

   end if;


end Tanh;