separate ( Generic_Elementary_Functions )

function Exp( X : Float_Type ) return Float_Type is

-- On input, X is a floating-point value in Float_Type;
-- On output, the value of exp(X) (the natural exponential of X) is returned.

-- 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.

   Result                             : Float_Type;

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

   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 :=
      4.0 * Common_Float(Float_Type'Safe_Emax) * 6.931471806E-1;

   Small_Threshold  : constant Common_Float :=
      Float_Type'Base'Epsilon;
      

begin

-- Filter out exceptional cases.

   Y := Common_Float( X );

   if abs(Y) >= Large_Threshold then
      if ( Y > 0.0 ) then
         raise Constraint_Error;
      else
	 return( Float_Type( 0.0 ) );
      end if;
   elsif abs(Y) <= SMALL_THRESHOLD then
      return( Float_Type( 1.0 + Y ) );
   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 (Y, M, Z1, Z2);

   case Common_Float'Machine_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;

   Result := Float_Type( Scale( Y, M ) );

   return( Result );


end Exp;