separate ( Generic_Elementary_Functions )

function Arccosh( X : Float_Type ) return Float_Type is

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

-- The definition of Arccosh(Y) is log( Y + sqrt(Y*Y - 1) ), Y >= 1.
-- To obtain good accuracy, we consider several cases:
-- 1) Y = 1, simply return 0.
-- 2) 1 < Y <= sqrt(2)
--    Y + sqrt(Y*Y-1) = Y + sqrt( (Y-1)*(Y+1) )
-- 3) sqrt(2) < Y < 10/epsilon,
--    Y + sqrt(Y*Y-1) = 2( Y  -  0.5/[ sqrt(Y*Y-1) + Y ] ).
-- 4) 10/epsilon <= Y, then 
--    Y + sqrt(Y*Y-1) = 2Y for practical purposes.
-- Note that (3) and (4) are suited for invoking the kernel procedure
-- KP_Log(Input) which returns M, Z1, and Z2 where 
--    log(Input) = M * log(2)   +   Z1   +   Z2.
--
  
   Y1, Y, M, Z1, Z2, Result                   : Common_Float;

   Zero  : constant := 0.0;
   One   : constant := 1.0;
   Half  : constant := 0.5;
   Root2 : constant := 1.41421_35623_73095_04880_16887_24209_69807;
   
   Large_Threshold : constant Common_Float :=
		       10.0 / Common_Float'Base'Epsilon;
   Log2_Lead  : constant Common_Float := 16#0.B17#;
   Log2_Trail : constant Common_Float :=
               16#0.000217F7D1CF79ABC9E3B39803F2F6AF40#;

begin


   Y := Common_Float(X); 
   if (Y < One) then
      raise Argument_Error;
   end if;

   if (Y <= Root2) then

      if (Y = One) then
         return( Float_Type(Zero) );
      else
	 Y1 := Y - One;
         Y := Y1 + KF_Sqrt( Y1*(Y+One) );
         return( Float_Type(KF_L1p( Y )) );
      end if;

   else

      if (Y < Large_Threshold) then
         Y := Y - Half/(Y + KF_Sqrt( (Y-One)*(Y+One) ));
      end if;
      KP_Log( Y, M, Z1, Z2 );
      M := M + One;

   end if;

   Result := M*Log2_Lead + (Z1 + (Z2 + M*Log2_Trail));


   return( Float_Type(Result) );

end Arccosh;