|
DataMuseum.dkPresents historical artifacts from the history of: Rational R1000/400 Tapes |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about Rational R1000/400 Tapes Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - downloadIndex: ┃ B T ┃
Length: 12840 (0x3228)
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⟧
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;