separate ( Generic_Elementary_Functions )

function Arcsinh( X : Float_Type ) return Float_Type is

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

-- The definition of Arcsinh(Y) is log( Y + sqrt(Y*Y + 1) )
-- For symmetry, we return sign(Y)*Arcsinh(|Y|). The discussions below
-- therefore assume Y >= 0.
-- To obtain good accuracy, we consider several cases:
-- 1) Y <= epsilon, simply return Y.
-- 2) epsilon < Y <= 0.5, 
--    Y + sqrt(Y*Y+1) = 1 + ( Y  +  Y*Y/[1 + sqrt(Y*Y+1)] )
--                    = 1 + ( Y  +  Y / [ (1/Y) + sqrt(1 + [1/Y*Y]) ] )
--    A formula best suited for the kernel function L1p.
-- 3) 0.5 < 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.
--
  
   Y, Sign_X, V, M, Z1, Z2, Result                   : Common_Float;

   One  : constant := 1.0;
   Half : constant := 0.5;

   Small_Threshold : constant Common_Float := 
				 Common_Float'Base'Epsilon;
   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); Sign_X := Copy_Sign(One,Y); Y := abs(Y);

   if (Y <= Half) then

      if (Y < Small_Threshold) then
         return( Float_Type( Copy_Sign(Y,Sign_X) ) );
      else
         V := Y + Y/( (One/Y) + KF_Sqrt(One + (One/(Y*Y))) );
         Result := KF_L1p( V );
      end if;

   else

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

   end if;

   return( Float_Type( Copy_Sign( Result, Sign_X ) ) );

end Arcsinh;