separate ( Generic_Elementary_Functions )

function Arctanh( X : Float_Type ) return Float_Type is

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

-- The definition of Arctanh(Y) is log((1+Y)/(1-Y)) / 2, which is also 
-- equivalent to the following three formulas:
--      1.  ( log(1+Y) - log(1-Y) ) / 2
--      2.  ( log(Y+1) - log(-Y+1) ) / 2.
--      3.  log( 1 + ( (2*Y) / (1-Y) ) ) / 2.
-- but computationally, the last formula is better.

   Z, Sign_Y                            : Common_Float;

   Y, Abs_Y, Temp                       : Common_Float;

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

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


begin

-- Filter out exceptional cases.

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

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

   if ( Abs_Y = 1.0 ) then
      raise Constraint_Error;
   end if;

   if ( Abs_Y > 1.0 ) then
      raise Argument_Error;
   end if;

-- Calculate Arctanh(Y) by using KF_L1p.

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

   Temp   := ( 2.0 * Abs_Y ) / ( 1.0 - Abs_Y );
   Temp   := KF_L1p ( Temp );
   Z      := Sign_Y *  0.5 * Temp;
   return ( Float_Type(Z) );

end Arctanh;