separate ( Generic_Elementary_Functions )

procedure KP_Exp( Y      : in  Common_Float;
                  M      : out Common_Int;
                  Z1, Z2 : out Common_Float) is

-- On input, Y is a floating-point value in Common_Float;
-- On output, the values M, Z1, and Z2 are returned such that
--      exp(Y) = 2**M * (Z1 + Z2),  where
-- exp(Y) is the natural exponential of Y.

-- Computations needed in order to obtain M, Z1, and Z2 involve three steps.

-- First, we reduce the argument Y to the form
--          Y = N * log2/32 + remainder,
-- where N has the value of an integer and |remainder| <= log2/64.
-- The value of N = Y * 32/log2 rounded to the nearest integer and
-- the remainder = Y - N*log2/32.

-- Second, we approximate exp(Y) - 1 by using KF_Em1Sm (Y1, Y2)
-- where Y1 is the leading part of the remainder and Y2 is the trailing part
-- of the remainder.  The approximation is carried out by using a polynomial
-- and Working_Float.  See the code in KF_Em1Sm for more details.

-- Third, we reconstruct the exponential of Y so that
--   exp(Y) = 2**M * (Z1 + Z2).

   F, X, R1, R2, Q, Tmp        : Common_Float;

   Two_to_6            : constant := 64;
   Two_to_8            : constant := 256;

   Log2_by_32          : constant :=
      16#5.8B90B_FBE8E_7BCD5_E4F1D_9CC01_F97B5_7A079_A1933_94C#E-2;

   Thirty_Two_by_Log2  : constant :=
      16#2E.2A8EC_A5705_FC2EE_FA1FF_B41A4_74FA2_3AD5E#;

   N, N_1, N_2, J              : Common_Int;

begin

   X := Y;

-- Step 1. Reduce the argument.

-- To perform argument reduction, we find the integer N such that
--   X = N * log2/32 + remainder, |remainder| <= log2/64.
-- N is defined by round-to-nearest-integer( X*32/log2 ) and
-- remainder by X - N*log2/32. The calculation of N is
-- straightforward whereas the computation of X - N*log2/32
-- must be carried out carefully.
-- log2/32 is so represented in two pieces that
-- (1) log2/32 is known to extra precision, (2) the product
-- of N and the leading piece is a model number and is hence
-- calculated without error, and (3) the subtraction of the value
-- obtained in (2) from X is a model number and is hence again obtained
-- without error.

-- Perform argument reduction in Common_Float.

   N := Common_Int( X * Thirty_Two_by_Log2 );
   if abs(N) >= Two_to_8 then
      N_2 := N mod Two_to_6;
      N_1 := N - N_2;
   else
      N_2 := 0;
      N_1 := N;
   end if;
   Tmp := Common_Float( N_1 ) * Log2_by_32_Lead;
   if abs( X ) >= abs( Tmp ) then
      R1 := X - Tmp;
   else
      Tmp := 0.5 * Tmp;
      R1 := (X - Tmp) - Tmp;
   end if;
   if N_2 /= 0 then
      R1 := R1 - Common_Float(N_2) * Log2_by_32_Lead;
   end if;
   R2 := -Common_Float(N) * Log2_by_32_Trail;
   J := N mod 32;
   M := (N - J)/32;

-- Step 2. Approximation.

   Q := KF_Em1Sm( R1, R2 );

-- Step 3. Function value reconstruction.

-- We now reconstruct the exponential of the input argument
-- so that Exp(Y) = 2**M * (Z1 + Z2).
-- The order of the computation below must be strictly observed.

   F  := Two_to_Jby32( J, Lead ) + Two_to_Jby32( J, Trail );

   Z1 := Two_to_Jby32( J, Lead );
   Z2 := Two_to_Jby32( J, Trail ) + F*Q;

   return;


end KP_Exp;