DataMuseum.dk

Presents historical artifacts from the history of:

Rational R1000/400 Tapes

This is an automatic "excavation" of a thematic subset of
artifacts from Datamuseum.dk's BitArchive.

See our Wiki for more about Rational R1000/400 Tapes

Excavated with: AutoArchaeologist - Free & Open Source Software.


top - download
Index: ┃ B T

⟦c6bcf374e⟧ TextFile

    Length: 12840 (0x3228)
    Types: TextFile
    Names: »B«

Derivation

└─⟦5f3412b64⟧ Bits:30000745 8mm tape, Rational 1000, ENVIRONMENT 12_6_5 TOOLS 
    └─ ⟦91c658230⟧ »DATA« 
        └─⟦458657fb6⟧ 
            └─⟦1472c4407⟧ 
                └─⟦this⟧ 
└─⟦d10a02448⟧ Bits:30000409 8mm tape, Rational 1000, ENVIRONMENT, D_12_7_3
    └─ ⟦fc9b38f02⟧ »DATA« 
        └─⟦9b46a407a⟧ 
            └─⟦2e03b931c⟧ 
                └─⟦this⟧ 

TextFile

separate (Generic_Elementary_Functions)

procedure Kp_Logpw (X : in Common_Float; Z1, Z2 : out Common_Float) is

-- On input, X is a floating-point value in Common_Float;
-- On output, the values Z1, and Z2 are returned such that
--      log(X) =  Z1 + Z2 , where
-- log(X) is the natural log of X. Moreover, Z1 only has a few
-- digits followed by trailing zeros.

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

-- First, we decompose the argument X to the form
--          X  =  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.
--
-- It is not hard to see that
--    log(X) = M*log(2) + log(F1) + log( 1 + F2/F1 ).
-- Hence, we return Z1 + Z2 as an accurate approximation to the
-- sum above. The idea is that Z1 is the first few digits of
-- the sum. Z2 is the correction to the sum accurate to full machine
-- precision.
--
--

-- If X is close to 1, satisfying |log(X)| < 1/64, then we use
-- an approximation similar to the core approximation in Step 2
-- above to calculate log( 1 + (X-1) ).

-- 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.

   Y, F, F1, F2, G, U, U1, U2, V, W1, W2, Half_F2,
   Half_F2_H, Half_F2_L, Poly, A1, A2, T1, T2 : Common_Float;
   Index, I, K : Common_Int;
   L, Ll : Positive;
   Thres1 : constant := 16#0.FC07F560#;
   Thres2 : constant := 16#1.04080AB5#;

begin

   L  := (2 * Common_Float'Machine_Mantissa) / 5;
   Ll := Common_Float'Machine_Mantissa;

-- I. Filter out exceptional arguments.

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


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

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

      if (X >= 1.0) then
         F := X - 1.0;
      else
         F := X - 0.5;
         F := F - 0.5;
      end if;
-- On machines with guard digit for subtraction, the conditions
-- above can be replaced by   X := X - 1.0


      F1 := Leading_Part (F, L);
      F2 := F - F1;
      G  := 1.0 / (2.0 + F);
      U  := 2.0 * F * G;
      V  := U * U;
      U1 := Leading_Part (U, L);
      U2 := G * ((2.0 * (F - U1) - U1 * F1) - U1 * F2);

-- We now compute a polynomial approximation to log(1 + F2/F1)
      case Float_Type'Base'Digits is
         when 1 .. 6 =>
            declare
               type Working_Float is digits 6;
               Q, R, S : Working_Float;
            begin
               R  := Working_Float (U);
               S  := Working_Float (V);
               Q  := S * 8.33340_08285_51364E-02;
               Z2 := U2 + Common_Float (R * Q);
            end;

         when 7 .. 15 =>
            declare
               type Working_Float is digits (15 + System.Max_Digits -
                                             abs (15 - System.Max_Digits)) / 2;
               -- this is min( 15, System.Max_Digits )
               Q, R, S : Working_Float;
            begin
               R := Working_Float (U);
               S := Working_Float (V);
               Q := S * (8.33333_33333_33335_93622E-02 +
                         S * (1.24999_99997_81386_68903E-02 +
                              S * (2.23219_81075_85598_51206E-03)));

               Z2 := U2 + Common_Float (R * Q);
            end;

         when 16 =>
            declare
               type Working_Float is digits (16 + System.Max_Digits -
                                             abs (16 - System.Max_Digits)) / 2;
               Q, R, S : Working_Float;
            begin
               R := Working_Float (U);
               S := Working_Float (V);
               Q := S * (8.33333_33333_33335_93622E-02 +
                         S * (1.24999_99997_81386_68903E-02 +
                              S * (2.23219_81075_85598_51206E-03)));

               Z2 := U2 + Common_Float (R * Q);
            end;

         when 17 .. 18 =>
            declare
               type Working_Float is digits (18 + System.Max_Digits -
                                             abs (18 - System.Max_Digits)) / 2;
               Q, R, S : Working_Float;
            begin
               R := Working_Float (U);
               S := Working_Float (V);
               Q := S * (8.33333_33333_33335_93622E-02 +
                         S * (1.24999_99997_81386_68903E-02 +
                              S * (2.23219_81075_85598_51206E-03)));

               Z2 := U2 + Common_Float (R * Q);
            end;

         when 19 .. 27 =>
            declare
               type Working_Float is digits (27 + System.Max_Digits -
                                             abs (27 - System.Max_Digits)) / 2;
               Q, R, S : Working_Float;
            begin
               R  := Working_Float (U);
               S  := Working_Float (V);
               Q  :=
                  S *
                     (8.33333_33333_33333_33333_33334_07301_529E-02 +
                      S *
                         (1.24999_99999_99999_99998_61732_74718_869E-02 +
                          S *
                             (2.23214_28571_42866_13712_34336_23012_985E-03 +
                              S *
                                 (4.34027_77751_26439_67391_35491_00214_979E-04 +
                                  S *
                                     (8.87820_39767_24501_02052_39367_49695_054E-05)))));
               Z2 := U2 + Common_Float (R * Q);
            end;

         when 28 .. 33 =>
            declare
               type Working_Float is digits (33 + System.Max_Digits -
                                             abs (33 - System.Max_Digits)) / 2;
               Q, R, S : Working_Float;
            begin
               R  := Working_Float (U);
               S  := Working_Float (V);
               Q  :=
                  S *
                     (8.33333_33333_33333_33333_33333_33332_96298_39318E-02 +
                      S *
                         (1.25000_00000_00000_00000_00000_93488_19499_40702E-02 +
                          S *
                             (2.23214_28571_42857_14277_26598_59261_40273_30694E-03 +
                              S *
                                 (4.34027_77777_77814_30973_20354_95180_362E-04 +
                                  S *
                                     (8.87784_09009_03777_78533_78449_15942_610E-05 +
                                      S *
                                         (1.87809_65740_24066_11924_19609_24471_232E-05))))));
               Z2 := U2 + Common_Float (R * Q);
            end;

         when others =>
            raise Program_Error;  -- assumption (1) is violated.

      end case;

      Z1 := U1;
      return;

   end if;
-- End of the case when X is close to 1

-- 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, K);

   if (Radix = 16) then
      K := 4 * K;
      while F < 0.5 loop
         F := F + F;
         K := K - 1;
      end loop;
   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;
-- 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 )
   G         := 1.0 / (F1 + Half_F2);
   U         := F2 * G;
   V         := U * U;
   U1        := Leading_Part (U, L);
   Half_F2_H := Leading_Part (Half_F2, L);
   Half_F2_L := Half_F2 - Half_F2_H;
   U2        := G * (((F2 - U1 * F1) - U1 * Half_F2_H) - U1 * Half_F2_L);

-- Approximate [log(1+U/2)-log(1-U/2)]/U  -   1.
   case Float_Type'Base'Digits is
      when 1 .. 6 =>
         declare
            type Working_Float is digits 6;
            Q, R, S : Working_Float;
         begin
            R  := Working_Float (U);
            S  := Working_Float (V);
            Q  := S * 8.33340_08285_51364E-02;
            W2 := U2 + Common_Float (R * Q);
         end;

      when 7 .. 15 =>
         declare
            type Working_Float is
               digits (15 + System.Max_Digits - abs (15 - System.Max_Digits)) /
                      2;
            -- this is min( 15, System.Max_Digits )
            Q, R, S : Working_Float;
         begin
            R := Working_Float (U);
            S := Working_Float (V);
            Q := S * (8.33333_33333_33335_93622E-02 +
                      S * (1.24999_99997_81386_68903E-02 +
                           S * (2.23219_81075_85598_51206E-03)));

            W2 := U2 + Common_Float (R * Q);
         end;

      when 16 =>
         declare
            type Working_Float is
               digits (16 + System.Max_Digits - abs (16 - System.Max_Digits)) /
                      2;
            Q, R, S : Working_Float;
         begin
            R := Working_Float (U);
            S := Working_Float (V);
            Q := S * (8.33333_33333_33335_93622E-02 +
                      S * (1.24999_99997_81386_68903E-02 +
                           S * (2.23219_81075_85598_51206E-03)));

            W2 := U2 + Common_Float (R * Q);
         end;

      when 17 .. 18 =>
         declare
            type Working_Float is
               digits (18 + System.Max_Digits - abs (18 - System.Max_Digits)) /
                      2;
            Q, R, S : Working_Float;
         begin
            R := Working_Float (U);
            S := Working_Float (V);
            Q := S * (8.33333_33333_33335_93622E-02 +
                      S * (1.24999_99997_81386_68903E-02 +
                           S * (2.23219_81075_85598_51206E-03)));

            W2 := U2 + Common_Float (R * Q);
         end;

      when 19 .. 27 =>
         declare
            type Working_Float is
               digits (27 + System.Max_Digits - abs (27 - System.Max_Digits)) /
                      2;
            Q, R, S : Working_Float;
         begin
            R  := Working_Float (U);
            S  := Working_Float (V);
            Q  :=
               S *
                  (8.33333_33333_33333_33333_33334_07301_529E-02 +
                   S *
                      (1.24999_99999_99999_99998_61732_74718_869E-02 +
                       S *
                          (2.23214_28571_42866_13712_34336_23012_985E-03 +
                           S *
                              (4.34027_77751_26439_67391_35491_00214_979E-04 +
                               S *
                                  (8.87820_39767_24501_02052_39367_49695_054E-05)))));
            W2 := U2 + Common_Float (R * Q);
         end;

      when 28 .. 33 =>
         declare
            type Working_Float is
               digits (33 + System.Max_Digits - abs (33 - System.Max_Digits)) /
                      2;
            Q, R, S : Working_Float;
         begin
            R  := Working_Float (U);
            S  := Working_Float (V);
            Q  :=
               S *
                  (8.33333_33333_33333_33333_33333_33332_96298_39318E-02 +
                   S *
                      (1.25000_00000_00000_00000_00000_93488_19499_40702E-02 +
                       S *
                          (2.23214_28571_42857_14277_26598_59261_40273_30694E-03 +
                           S *
                              (4.34027_77777_77814_30973_20354_95180_362E-04 +
                               S *
                                  (8.87784_09009_03777_78533_78449_15942_610E-05 +
                                   S *
                                      (1.87809_65740_24066_11924_19609_24471_232E-05))))));
            W2 := U2 + Common_Float (R * Q);
         end;

      when others =>
         raise Program_Error;  -- assumption (1) is violated.

   end case;

   A2 := Common_Float (K);
   A1 := A2 * Pow_Log2_Lead + Pow_Log_Of_Jby64 (Index, Lead);
   A2 := A2 * Pow_Log2_Trail + Pow_Log_Of_Jby64 (Index, Trail);
   T1 := Leading_Part (A1 + U1, Ll - 2);     -- force storage
   T2 := (A1 - T1) + U1;
   T2 := T2 + A2;
   W1 := Leading_Part (T1 + W2, L);
   Z2 := ((T1 - W1) + W2) + T2;
   Z1 := W1;

   return;

end Kp_Logpw;