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.

   F, F1, F2, G, U, U1, U2, V, W1, W2,
   Half_F2, Half_F2_H, Half_F2_L,
   A1, A2, T1, T2                       : Common_Float;
   Index, 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;