separate ( Generic_Elementary_Functions ) function Cosh( X : Float_Type ) return Float_Type is -- On input, X is a floating-point value in Float_Type; -- On output, the value of cosh(X) (the hyperbolic cosine of X) is returned. -- The definition of cosh(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 : Common_Float; Y, Abs_Y, Z1, Z2 : Common_Float; M, J : Common_Int; 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; Cond : constant Common_Float := Base_Digits * Log2; begin -- Filter out exceptional cases. Y := Common_Float( X ); Abs_Y := abs(Y); if Abs_Y >= Large_Threshold then raise Constraint_Error; end if; -- 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; -- Now, Z = 1/2 * exp( abs(X) ). Z := Scale( Y, M ); if (Abs_Y >= Cond) then -- When abs(Y) gets so big, adding (1/4)/Z will not make a difference in the -- outcome of cosh(X). return ( Float_Type(Z) ); else return ( Float_Type( Z + 0.25/Z ) ); end if; end Cosh;