⟦942430b8e⟧

E3 Source Code



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#;

   I, 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;
DataMuseum.dk

Rational R1000/400

⟦942430b8e⟧ Ada Source

Derivation

E3 Source Code

E3 Meta Data
