|
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 - download
Length: 9216 (0x2400)
Types: Ada Source
Notes: 03_class, FILE, R1k_Segment, e3_tag, function Kf_R1psm, seg_0130d6, separate Generic_Elementary_Functions
└─⟦8527c1e9b⟧ Bits:30000544 8mm tape, Rational 1000, Arrival backup of disks in PAM's R1000
└─ ⟦5a81ac88f⟧ »Space Info Vol 1«
└─⟦this⟧
separate (Generic_Elementary_Functions)
function Kf_R1psm (Y : Common_Float) return Common_Float is
-- On input, Y lies in the interval [-1/64, 1/64]
-- On output, the value sqrt( 1 + Y ) - 1 is returned.
-- The approximation employs both minmax polynomial and Newton Iteration.
--
-- First, sqrt(1+Y)-1 is approximated by a polynomial of the form
-- 0.5 Y - Y^2*(A1 - Y*(A2 - Y*(A3 ... ))).
-- Depending of the number of digits required, this polynomial may or may
-- not be the value return. Should the polynomial fall short of the
-- accuracy requirement, one step of the Newton iteration is applied to
-- the equation
-- (P + 1)^2 = 1 + Y, where P is the unknown.
-- Thus,
-- new P := P - (P^2 + 2P - Y)/(2P + 2)
-- where the initial guess P is the value of the polynomial.
Poly : Common_Float;
begin
case Float_Type'Base'Digits is
when 1 .. 6 =>
declare
type Working_Float is digits 6;
R, P : Working_Float;
begin
R := Working_Float (Y);
P := R * 0.5 - R * R * (0.12500_79870 - R * (0.62505_22999E-01));
Poly := Common_Float (P);
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 )
R, P, Two_P, Correction : Working_Float;
begin
R := Working_Float (Y);
P := R * 0.5 - R * R * (0.12500_79870 - R * (0.62505_22999E-01));
Two_P := P + P;
Correction := (P * P + (Two_P - R)) / (Two_P + 2.0);
Poly := Common_Float (P) - Common_Float (Correction);
end;
when 16 =>
declare
type Working_Float is
digits (16 + System.Max_Digits - abs (16 - System.Max_Digits)) /
2;
R, P, Two_P, Correction : Working_Float;
begin
R := Working_Float (Y);
P := R * 0.5 - R * R * (0.12500_79870 - R * (0.62505_22999E-01));
Two_P := P + P;
Correction := (P * P + (Two_P - R)) / (Two_P + 2.0);
Poly := Common_Float (P) - Common_Float (Correction);
end;
when 17 .. 18 =>
declare
type Working_Float is
digits (18 + System.Max_Digits - abs (18 - System.Max_Digits)) /
2;
R, P, Two_P, Correction : Working_Float;
begin
R := Working_Float (Y);
P := R * 0.5 -
R * R * (0.125 - R * (0.62505_85095_32379_40606E-01 -
R * (0.39067_55576_62793_58660E-01)));
Two_P := P + P;
Correction := (P * P + (Two_P - R)) / (Two_P + 2.0);
Poly := Common_Float (P) - Common_Float (Correction);
end;
when 19 .. 27 =>
declare
type Working_Float is
digits (27 + System.Max_Digits - abs (27 - System.Max_Digits)) /
2;
R, P, Two_P, Correction : Working_Float;
begin
R := Working_Float (Y);
P :=
R * 0.5 -
R * R *
(0.12500_00000_00043_29869_79603 -
R *
(0.62499_99949_58313_28657_64839E-01 -
R *
(0.39062_49897_81141_25807_82311E-01 -
R *
(0.27349_66395_87080_51587_13279E-01 -
R * (0.20514_47322_06561_78273_82000E-01)))));
Two_P := P + P;
Correction := (P * P + (Two_P - R)) / (Two_P + 2.0);
Poly := Common_Float (P) - Common_Float (Correction);
end;
when 28 .. 33 =>
declare
type Working_Float is
digits (33 + System.Max_Digits - abs (33 - System.Max_Digits)) /
2;
R, P, Two_P, Correction : Working_Float;
begin
R := Working_Float (Y);
P :=
R * 0.5 -
R * R *
(0.12500_00000_00035_79677_13859_89202_617 -
R *
(0.62500_00000_00664_48622_98045_60336E-01 -
R *
(0.39062_49911_61750_24807_32583_41508E-01 -
R *
(0.27343_74885_03551_93697_91233_69204E-01 -
R *
(0.20514_00699_96599_52762_98109_70261E-01 -
R *
(0.16119_54134_09746_99174_73737_89068E-01))))));
Two_P := P + P;
Correction := (P * P + (Two_P - R)) / (Two_P + 2.0);
Poly := Common_Float (P) - Common_Float (Correction);
end;
when others =>
raise Program_Error; -- assumption (1) is violated.
end case;
return (Poly);
end Kf_R1psm;
nblk1=8
nid=0
hdr6=10
[0x00] rec0=1f rec1=00 rec2=01 rec3=056
[0x01] rec0=1d rec1=00 rec2=02 rec3=006
[0x02] rec0=00 rec1=00 rec2=08 rec3=00c
[0x03] rec0=18 rec1=00 rec2=03 rec3=04c
[0x04] rec0=19 rec1=00 rec2=04 rec3=032
[0x05] rec0=00 rec1=00 rec2=07 rec3=03a
[0x06] rec0=19 rec1=00 rec2=05 rec3=030
[0x07] rec0=12 rec1=00 rec2=06 rec3=001
tail 0x2150db5be82b151add51b 0x42a00066462061e03