DataMuseum.dk

top - download

    Length: 5386 (0x150a)
    Types: TextFile
    Names: »B«

└─⟦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)

function Kf_L1psm (Y : Common_Float) return Common_Float is

-- On input, Y lies in the interval ( exp(-1/16)-1, exp(1/16)-1 ).
-- On output, the value log( 1 + Y ) is returned.
-- The approximation exploit the identity
--     log( 1 + Y ) = log( 1 + u/2 )  /  log( 1 - u/2 ), where
--     u  = 2Y / (2+Y).
-- Note that the right hand side has an odd Taylor series expansion
-- which converges much faster than the Taylor series expansion of
-- log( 1 + Y ) in Y. Thus, we approximate log( 1 + Y ) by
--     u + A1 * u^3 + A2 * u^5 + ... + An * u^(2n+1).

-- One subtlety is that since   u    cannot be calculated from
-- Y exactly, the rounding error in the first   u   should be
-- avoided if possible. To accomplish this, we observe that
--               u  =  Y  -  Y*Y/(2+Y).
-- Since Y is the input argument, and thus presumed exact, the
-- formula above approximates u accurately because
--               u  =  Y  -  correction,
-- and the magnitude of "correction" (of the order of Y*Y)
-- is small.
-- With these observations, we will approximate log( 1 + Y ) by
--    Y + (  (A1*u^3 + ... + An*u^(2n+1)) - correction ).

-- All calculations will be carried out in Common_Float.
-- Since this functions does not need to perform any argument
-- reduction or reconstruction, it will be faster than the
-- calculation of  log   in the normal case, despite computations
-- done solely in Common_Float.

   U, V, Correction, Poly : Common_Float;

begin

-- We approximate log(1+Y) by a odd polynomial in U, where
--          U = 2Y/(2+Y) = Y - Y*Y/(2+Y).
-- We now calculate U and other related quantities.

   U          := Y / (2.0 + Y);
   Correction := Y * U;
   U          := U + U;
   V          := U * U;

-- We use an appropriate polynomial approximation that will deliver
-- a result accurate enough with respect to Float_Type'Base'Digits.
-- Note that the upper bounds of the cases below (6, 15, 16, 18,
-- 27, and 33) are attributes of predefined floating types of
-- common systems.

   case Float_Type'Base'Digits is

      when 1 .. 6 =>

         Poly := U * V * (8.33333_18167E-02 + V * (1.25123_45996E-02));


      when 7 .. 15 =>

         Poly := U * V * (8.33333_33333_33179_23934E-02 +
                          V * (1.25000_00003_77175_09602E-02 +
                               V * (2.23213_99879_19448_06202E-03 +
                                    V * (4.34887_77770_76145_52256E-04))));

      when 16 =>

         Poly := U * V *
                    (8.33333_33333_33333_37909E-02 +
                     V * (1.24999_99999_99836_60424E-02 +
                          V * (2.23214_28590_37388_98829E-03 +
                               V * (4.34026_81922_36633_88536E-04 +
                                    V * (8.89985_05326_15867_46686E-05)))));

      when 17 .. 18 =>

         Poly :=
            U * V *
               (8.33333_33333_33333_37909_47954E-02 +
                V * (1.24999_99999_99836_60424_94790E-02 +
                     V * (2.23214_28590_37388_98829_86093E-03 +
                          V * (4.34026_81922_36633_88536_93268E-04 +
                               V * (8.89985_05326_15867_46686_13783E-05)))));

      when 19 .. 27 =>

         Poly :=
            U * V *
               (8.33333_33333_33333_33333_33333_32289_687E-02 +
                V *
                   (1.25000_00000_00000_00000_00085_22215_240E-02 +
                    V *
                       (2.23214_28571_42857_14261_92866_48258_201E-03 +
                        V *
                           (4.34027_77777_77780_98809_33315_23070_295E-04 +
                            V *
                               (8.87784_09090_66940_14314_76792_89545_412E-05 +
                                V *
                                   (1.87800_48181_19495_76549_39975_36152_064E-05 +
                                    V *
                                       (4.06898_41092_88362_34597_46674_51012_081E-06 +
                                        V *
                                           (9.01142_71591_15642_85666_93997_67835_754E-07))))))));

      when 28 .. 33 =>

         Poly :=
            U * V *
               (8.33333_33333_33333_33333_33333_33333_62086_64017E-02 +
                V *
                   (1.24999_99999_99999_99999_99999_97094_37287_04079E-02 +
                    V *
                       (2.23214_28571_42857_14285_72439_65284_28894_60352E-03 +
                        V *
                           (4.34027_77777_77777_77605_29930_95879_387E-04 +
                            V *
                               (8.87784_09090_90925_73327_23371_16332_464E-05 +
                                V *
                                   (1.87800_48076_82610_80293_27069_90675_944E-05 +
                                    V *
                                       (4.06901_04514_65706_22856_93057_91655_559E-06 +
                                        V *
                                           (8.97568_30488_20751_96312_20868_52977_742E-07 +
                                            V *
                                               (2.01671_52613_96766_05830_17681_29488_302E-07)))))))));
      when others =>

         raise Program_Error;  -- assumption (1) is violated.

   end case;

-- return the result

   return (Y + (Poly - Correction));

end Kf_L1psm;