|
DataMuseum.dkPresents historical artifacts from the history of: Rational R1000/400 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about Rational R1000/400 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 8192 (0x2000)
Types: Ada Source
Notes: 03_class, FILE, R1k_Segment, e3_tag, function Kf_L1psm, seg_0130d3, separate Generic_Elementary_Functions
└─⟦8527c1e9b⟧ Bits:30000544 8mm tape, Rational 1000, Arrival backup of disks in PAM's R1000
└─⟦cfc2e13cd⟧ »Space Info Vol 2«
└─⟦this⟧
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;
nblk1=7
nid=0
hdr6=e
[0x00] rec0=16 rec1=00 rec2=01 rec3=04c
[0x01] rec0=1b rec1=00 rec2=02 rec3=056
[0x02] rec0=01 rec1=00 rec2=07 rec3=006
[0x03] rec0=1e rec1=00 rec2=03 rec3=044
[0x04] rec0=15 rec1=00 rec2=04 rec3=04e
[0x05] rec0=15 rec1=00 rec2=05 rec3=048
[0x06] rec0=12 rec1=00 rec2=06 rec3=000
tail 0x2170e744282b151a39bc8 0x42a00066462061e03