⟦c54f3bae8⟧

E3 Source Code



separate (Generic_Elementary_Functions)

procedure Kp_Log (Y : in Common_Float; M, 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
--      log(Y) = M*log(2) + Z1 + Z2 , where
-- log(Y) is the natural log of Y.

-- Except for arguments Y that are very close to 1, computations
-- needed in order to obtain M, Z1, and Z2 involve two steps.

-- First, we decompose the argument Y to the form
--          Y  =  2**M  *  (F1  +  F2),
-- where  1 <= F1+F2 < 2, M has the value of an integer,
-- F1 = 1 + j/64, j ranges from 0 to 64, and |F2| <= 1/128.

-- Second, we approximate log( 1 + F2/F1 ) by an odd polynomial
-- in U, where
--              U  =  2 F2 / (2 F2 + F1).
-- Note that
--     log( 1 + F2/F1 ) = log( 1 + U/2 ) - log( 1 - U/2 ).
-- The core approximation calculates
--     Poly = [log( 1 + U/2 ) - log( 1 - U/2 )]/U   -   1.
-- Note that  log(1 + U/2) - log(1 - U/2) = 2 arctanh ( U/2 ),
-- thus,
--     Poly =  2 arctanh( U/2 ) / U    -     1.
-- That value is computed by the kernel function   KF_Atanh.

-- It is not hard to see that
--    log(Y) = M*log(2) + log(F1) + log( 1 + F2/F1 ).
-- Hence, we return Z1 := log(F1), and  Z2 := log( 1 + F2/F1).
-- The values of log(F1) are calculated beforehand and stored in the program.

-- If Y is close to 1, satisfying |log(Y)| < 1/16, then we use
-- an approximation similar to the core approximation in Step 2
-- above to calculate log( 1 + (Y-1) ). For details of the
-- computation, refer to the code and comments in the separate
-- body of the function KF_L1pSm.

-- As far as the types used for internal computation are concerned,
-- we use Common_Float to carry out Step 1,  and use
-- Working_Float to carry out Step 2. See the code in the separate
-- body for the function KF_Atanh for details about Working_Float.

   X, F, F1, F2, U, K, Half_F2, Poly : Common_Float;
   Index, I : Common_Int;
   Thres1 : constant := 16#0.F07D5#;
   Thres2 : constant := 16#1.1082C#;

   Log2_Lead  : constant Common_Float := 16#0.B17#;
   Log2_Trail : constant Common_Float :=
      16#0.000217F7D1CF79ABC9E3B39803F2F6AF40#;

begin


-- I. Filter out exceptional arguments.

   if (Y <= 0.0) then
      raise Argument_Error;
   end if;

   X := Y;

-- II. Check if the input argument is very close to 1.0

   if (X >= Thres1 and X <= Thres2) then

-- Invoke the function KF_L1pSm to compute log(1 + (X-1))
-- for the value of Z1 and set M and Z2 to zero.
-- Although the value of X-1 is representable, on machines
-- without a guard digit for subtraction, X-1 may not be calculated
-- exactly.  Consequently, we have to take extra precaution.
      if (X >= 1.0) then
         X := X - 1.0;
      else
         X := X - 0.5;
         X := X - 0.5;
      end if;
-- On machines with guard digit for subtraction, the conditions
-- aove can be replaced by   X := X - 1.0

      M  := 0.0;
      Z1 := Kf_L1psm (X);
      Z2 := 0.0;
      return;

   end if;

-- III. General case.

-- Step 1. Decomposition.

-- Use the primitive function to extract K and F where
--     X  =  radix**K   *   F,   1/radix <= F < 1.
   Decompose (X, F, I);
   K := Common_Float (I);

   if (Common_Float'Machine_Radix = 16) then
      K := 4.0 * K;
      while F < 0.5 loop
         F := F + F;
         K := K - 1.0;
      end loop;
   elsif (Common_Float'Machine_Radix /= 2) then
      raise Program_Error; -- Assumption (1) is violated.
   end if;
-- Now  X = 2**K  * F,  1/2 <= F < 1.

   Index := Common_Int (2#1.0#E7 * F);
   F1    := Common_Float (Index) * 2#1.0#E-7;
-- Again, on machines with gaurd digit for subtraction, the
-- statements below can be replaced by  Half_F2 := F - F1.
   if (F >= F1) then
      Half_F2 := F - F1;
   else
      F2      := F1 * 0.5;
      Half_F2 := (F - F2) - F2;
   end if;

   F1 := F1 + F1;
   F2 := Half_F2 + Half_F2;
   K  := K - 1.0;
-- At this point, X = 2**K * ( F1  +  F2 ) where
-- F1 := j/64, j = 64, 65, ..., 128 and |F2| <= 1/128.

-- Step 2. Approximation.

-- Calculate U := 2 F2 / ( 2 F1 + F2 ) = F2 / ( F1 + Half_F2 )
   U := F2 / (F1 + Half_F2);

-- Approximate [log(1+U/2)-log(1-U/2)]/U  -   1.
   Poly := Kf_Atanh (U);

   M  := K;
   Z1 := Log_Of_Jby64 (Index, Lead);
   Z2 := U + (Log_Of_Jby64 (Index, Trail) + U * Poly);
   return;

end Kp_Log;
DataMuseum.dk

Rational R1000/400

⟦c54f3bae8⟧ Ada Source

Derivation

E3 Source Code

E3 Meta Data
