⟦605182f8a⟧

E3 Source Code



separate (Generic_Elementary_Functions)

function Kf_Logb (X, Base : Float_Type) return Float_Type is

-- On input, X and Base are floating-point values in Float_Type;
-- on output, the value of log(X) (log of X to the Base) is returned.

-- The bulk of the computations are performed by the procedure
-- KP_Log( Y, M, Z1, Z2 ) which returns log(Y) in M, Z1, and Z2
-- where
--               log(Y) = M * log2  +  Z1  +  Z2,
-- M of integer value, and Z1 only has at most 12 significant bits.

-- Two special cases are considered : Base = 2.0 and Base = 10.0.
-- Since  log_2(X) = log(X) * (1/log(2)) and
--        log(X)   = M log 2 + Z1 + Z2.
-- Thus,  log_2(X) = M  +  (Z1 + Z2) * (1/log(2)).
-- We store 1/log(2) as  Log2_Inv, Log2_Inv_Lead, and Log2_Inv_Trail.
-- If M = 0, we recast (Z1,Z2) as (W1,W2) so that W1*Log2_Inv_Lead is
-- exact. If M is non-zero, we simply compute
--        M  +  (Z1+Z2)*Log2_Inv.

-- Next, for log_10(X), we store Log10_Inv, Log10_Inv_Lead, and
-- Log10_Inv_Trail. The trick is the same for M = 0. For M /= 0,
-- we need to compute  M * log(2)/log(10)  +  (Z1+Z2)/log(10).
-- The constants  Log2div10_Lead and Log2div10_Trail are the leading
-- and trailing part of log(2)/log(10).


   Y, Z, M, Z1, Z2, W1, W2, Temp, Result : Common_Float;
   Ndigits : Integer := Common_Float'Machine_Mantissa / 2;

   Log2_Lead  : constant Common_Float := 16#0.B17#;
   Log2_Trail : constant Common_Float :=
      16#0.00021_7F7D1_CF79A_BC9E3_B3980_3F2F6_AF40#;

   Log2_Inv_Lead  : constant Common_Float := 16#1.71#;
   Log2_Inv_Trail : constant Common_Float :=
      16#0.00547_652B8_2FE17_77D0F_FDA0D_23A7D_11D6A_EF55#;
   Log2_Inv       : constant Common_Float :=
      16#1.71547_652B8_2FE17_77D0F_FDA0D_23A7D_11D6A_EF55#;

   Log10_Inv_Lead  : constant Common_Float := 16#0.6F2#;
   Log10_Inv_Trail : constant Common_Float :=
      16#0.000DE_C549B_9438C_A9AAD_D557D_699EE_191F7_1A301#;
   Log10_Inv       : constant Common_Float :=
      16#0.6F2DE_C549B_9438C_A9AAD_D557D_699EE_191F7_1A301#;

   Log2div10_Lead  : constant Common_Float := 16#0.4D1#;
   Log2div10_Trail : constant Common_Float :=
      16#0.00004_D427D_E7FBC_C47C4_ACD60_5BE48_BC135_69862#;

   Log16_Inv_Lead  : constant Common_Float := 16#0.5C5#;
   Log16_Inv_Trail : constant Common_Float :=
      16#0.00051_D94AE_0BF85_DDF43_FF683_48E9F_4475A_BBD54#;
   Log16_Inv       : constant Common_Float :=
      16#0.5C551_D94AE_0BF85_DDF43_FF683_48E9F_4475A_BBD54#;


begin

   Y := Common_Float (X);

   if (Base <= 0.0 or Base = 1.0) then
      raise Argument_Error;
   end if;

   if (Y = 0.0) then
      raise Constraint_Error;
   end if;

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

-- Get values of M, Z1, and Z2 so that the log of X to the base
-- can be calculated by log (x) = log(x) / log(base)

   Kp_Log (Y, M, Z1, Z2);

   if (Base = 2.0) then

      W1 := Leading_Part (Z1 + Z2, Ndigits);
      W2 := (Z1 - W1) + Z2;
      W2 := W1 * Log2_Inv_Trail + W2 * Log2_Inv;
      W1 := W1 * Log2_Inv_Lead;
      if (M = 0.0) then
         Result := W1 + W2;
      else
         Z1     := M + W1;
         Z2     := ((M - Z1) + W1) + W2;
         Result := Z1 + Z2;
      end if;

   elsif (Base = 10.0) then

      W1 := Leading_Part (Z1 + Z2, Ndigits);
      W2 := (Z1 - W1) + Z2;
      W2 := (W1 * Log10_Inv_Trail + W2 * Log10_Inv);
      W1 := W1 * Log10_Inv_Lead;

      if (M = 0.0) then
         Result := W1 + W2;
      else
         Y      := M * Log2div10_Lead;
         Z1     := Y + W1;
         Z2     := ((Y - Z1) + W1) + (W2 + M * Log2div10_Trail);
         Result := Z1 + Z2;
      end if;

   elsif (Base = 16.0) then

      W1 := Leading_Part (Z1 + Z2, Ndigits);
      W2 := (Z1 - W1) + Z2;
      W2 := W1 * Log16_Inv_Trail + W2 * Log16_Inv;
      W1 := W1 * Log16_Inv_Lead;
      if (M = 0.0) then
         Result := W1 + W2;
      else
         M      := M * 0.25;
         Z1     := M + W1;
         Z2     := ((M - Z1) + W1) + W2;
         Result := Z1 + Z2;
      end if;

   else

      Z := Common_Float (Base);
      Kp_Log (Z, M, Z1, Z2);
      if (M = 0.0) then
         Temp := Z1 + Z2;
      else
         Temp := M * Log2_Trail + Z2;
         Temp := (M * Log2_Lead + Z1) + Temp;
      end if;
      Kp_Log (Y, M, Z1, Z2);
      if (M = 0.0) then
         Result := (Z1 + Z2) / Temp;
      else
         Result := M * Log2_Trail + Z2;
         Result := (M * Log2_Lead + Z1) + Result;
         Result := Result / Temp;
      end if;

   end if;
   return (Float_Type (Result));

end Kf_Logb;
DataMuseum.dk

Rational R1000/400

⟦605182f8a⟧ Ada Source

Derivation

E3 Source Code

E3 Meta Data
