separate ( Generic_Elementary_Functions )

function Sinh( X : Float_Type ) return Float_Type is

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

-- The definition of sinh(Y) is (exp(Y) - exp(-Y))/2, therefore
-- the bulk of the computations are performed by the procedure
-- KP_Exp (Y, M, Z1, Z2) which returns exp(Y) in M, Z1, and Z2
-- where
--              exp(Y) = 2**M * ( Z1 + Z2 )
-- M of integer value, and Z1 only has at most 12 significant bits.

   Z, Sign_Y                           : Common_Float;

   Y, Abs_Y, Z1, Z2, Cond              : Common_Float;
   M, J                                : Common_Int;

   One         : constant Common_Float := 1.0;

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

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

   Two_to      : constant array ( Common_Int range -3..3 ) of
                 Common_Float := ( 0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0 );

   Large_Threshold  : constant Common_Float :=
      8.0 * Common_Float(Float_Type'Safe_Emax) * 0.6931471806;


begin

-- Filter out exceptional cases.

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

   Y  := Common_Float( X );
   if ( Y > 0.0 ) then
      Sign_Y := One;
   else
      Sign_Y := -One;
   end if;

   Abs_Y  := abs(Y);


   if Abs_Y >= Large_Threshold then
      raise Constraint_Error;
   end if;

   Cond := Base_Digits * Log2;
   if (Abs_Y >= Cond) then
--    Get the values of M, Z1, and Z2  so that the natural exponential of Y
--    can be calculated by  Exp(Y) = 2**M * (Z1 + Z2)

      KP_Exp (Abs_Y, M, Z1, Z2);

      M := M - 1;

      case Radix is
         when 2 =>
	    Y      := Z1 + Z2;

         when others =>
	    J      := M rem 4;
	    M      := (M - J) / 4;
	    Z1     := Z1 * Two_to(J);
	    Z2     := Z2 * Two_to(J);
	    Y      := Z1 + Z2;
      end case;

      Z    := Sign_Y * Scale( Y, M );  

--    Now,  Z = sign(X) * 1/2 * exp( abs(X) ).
--    When abs(Y) gets so big, subtracting (1/4)/Z will not make
--    a difference in the outcome of the sinh(X).

      return ( Float_Type(Z) );

   else

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

   end if;

end Sinh;