⟦a1c086aef⟧

E3 Source Code



separate (Generic_Elementary_Functions)

function Kf_L1p (Y : Common_Float) return Common_Float is

-- On input, Y is a floating-point value in Common_Float;
-- On output, the value of log(1 + Y) is returned.

-- Except for arguments Y that are near zero, log( 1 + Y ) is
-- calculated by determining M, Z1, and Z2  such that
-- log( 1 + Y ) = M*Log2 + Z1 + Z2.  Then, except for arguments Y
-- that are very large numbers, computations needed in order to
-- retrieve the values of M, Z1, and Z2 involve two steps.

-- First, we decompose one plus the argument Y to the form
--          1 + 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 near zero, satisfying |log(Y)| < 1/16 - 1, then we use
-- an approximation similar to the core approximation in Step 2
-- above to calculate log( 1 + Y ). 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.

   F, F1, F2, U, M, Z, Z1, Z2, Half_F2, Poly, Result : Common_Float;
   Index, I, J : Common_Int;
   Thres1 : constant := 16#0.F07D5# - 1.0;
   Thres2 : constant := 16#1.1082C# - 1.0;

   Base_Digits : constant Common_Float :=
      Common_Float (3 * 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);

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

begin


-- I. Filter out exceptional case.

   if (Y <= -1.0) then
      raise Argument_Error;
   end if;

-- II. Check if the input argument is near zero.

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

-- Invoke the function KF_L1pSm to compute log(1 + Y).

      return (Kf_L1psm (Y));

   end if;

-- III. General case.

   if (Y >= (0.5 * Float_Type'Safe_Large)) then
     Kp_Log (Y, M, Z1, Z2);

   else

--    Calculate M, Z1, and Z2 not using KP_Log

--    Step 1.  Decomposition.

      Z := 1.0 + Y;

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

      if (Common_Float'Machine_Radix = 16) then

         M := 4.0 * M;
         while F < 0.5 loop
            J := J - 1;
            M := M - 1.0;
            F := F + F;
         end loop;
      end if;

--    Now, regardless if Radix is 2 or 16, we have
--       Z = 2**M * F = (radix**I) * (2**J) * F,  1/2 <= F < 1.

--    Extract the leading 7 bits of F.

      Index := Common_Int (F * 2#1.0#E7);
      F1    := Common_Float (Index) * 2#1.0#E-7;

--    We now calculate  F2  such that
--          2**M * (F1 + F2)                = 1 + Y  accurately.
--    i.e.  (radix**I) * (2**J) * (F1 + F2) = 1 + Y  accurately.

      if (M <= Base_Digits) then

--       F2 = 2**(-J) * ( ( r**(-I) - 2**(J) * F1 ) + r**(-I) * Y )

         if (J = 0) then
            F2 := Scale (1.0, -I) - F1;
            F2 := F2 + Scale (Y, -I);
         else
            F2 := Scale (1.0, -I) - Two_To (J) * F1;
            F2 := (F2 + Scale (Y, -I)) * Two_To (-J);
         end if;

      else

--       M is large.
--       F2 = 2**(-J) * ( ( r**(-I) * Y - (2**J) * F1 ) + r**(-I) )

         if (J = 0) then
            F2 := Scale (Y, -I) - F1;
            F2 := F2 + Scale (1.0, -I);
         else
            F2 := Scale (Y, -I) - Two_To (J) * F;
            F2 := (F2 + Scale (1.0, -I)) * Two_To (-J);
         end if;

      end if;

--    Now   2**M * (F1 + F2) = 1 + Y   accurately,
--    but   F1  = J/128, j = 64, ... , 128   and
--         |F2| < 1/256.

      M       := M - 1.0;
      F1      := F1 + F1;
      Half_F2 := F2;
      F2      := F2 + F2;

--    Now   2**M * (F1 + F2) = 1 + Y
--          F1  = j/64,  j = 64, ... , 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);

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

   end if;

-- Calculate log( 1 + Y ) from M, Z1, and Z2  such that
--    log(1 + Y ) = M*Log2 + Z1 + Z2.

   if (M = 0.0) then
      Result := Z1 + Z2;
   else
      Result := M * Log2_Trail + Z2;
      Result := (M * Log2_Lead + Z1) + Result;
   end if;

   return (Result);

end Kf_L1p;
DataMuseum.dk

Rational R1000/400

⟦a1c086aef⟧ Ada Source

Derivation

E3 Source Code

E3 Meta Data
