|
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 - metrics - downloadIndex: B T
Length: 16052 (0x3eb4)
Types: TextFile
Names: »B«
└─⟦a7d1ea751⟧ Bits:30000550 8mm tape, Rational 1000, !users!projects 94_04_11
└─⟦129cab021⟧ »DATA«
└─⟦9b477e385⟧
└─⟦f64eaa120⟧ Bits:30000752 8mm tape, Rational 1000, !projects 93 02 16
└─⟦6f12a12be⟧ »DATA«
└─⟦9b477e385⟧
└─⟦2f6cfab89⟧ Bits:30000547 8mm tape, Rational 1000, !projects 94-01-04
└─⟦d65440be7⟧ »DATA«
└─⟦9b477e385⟧
└─⟦this⟧
with Text_Io;
use Text_Io;
with Floating_Characteristics;
use Floating_Characteristics;
with Numeric_Io;
use Numeric_Io;
with Numeric_Primitives;
use Numeric_Primitives;
package body Core_Functions is
-- The following routines are coded directly from the algorithms and
-- coeficients given in "Software Manual for the Elementry Functions"
-- by William J. Cody, Jr. and William Waite, Prentice_Hall, 1980
-- CBRT by analogy
-- A more general formulation uses MANTISSA_TYPE, etc.
-- The coeficients are appropriate for 25 to 32 bits floating significance
-- They will work for less but slightly shorter versions are possible
-- The routines are coded to stand alone so they need not be compiled together
-- These routines have been coded to accept a general MANTISSA_TYPE
-- That is, they are designed to work with a manitssa either fixed of float
-- There are some explicit conversions which are required but these will
-- not cause any extra code to be generated
-- 16 JULY 1982 W A WHITAKER AFATL EGLIN AFB FL 32542
-- T C EICHOLTZ USAFA
function Sqrt (X : Float) return Float is
M, N : Exponent_Type;
F, Y : Mantissa_Type;
Result : Float;
subtype Index is Integer range 0 ..
100; -- #########################
Sqrt_L1 : Index := 3;
-- Could get away with SQRT_L1 := 2 for 28 bits
-- Using the better Cody-Waite coeficients overflows MANTISSA_TYPE
Sqrt_C1 : Mantissa_Type := 8#0.3317777777#;
Sqrt_C2 : Mantissa_Type := 8#0.4460000000#;
Sqrt_C3 : Mantissa_Type := 8#0.55202_36314_77747_36311_0#;
begin
if X = Zero then
Result := Zero;
return Result;
elsif X = One then
-- To get exact SQRT(1.0)
Result := One;
return Result;
elsif X < Zero then
New_Line;
Put ("CALLED SQRT FOR NEGATIVE ARGUMENT ");
raise Program_Error;
--Put (X);
Put (" USED ABSOLUTE VALUE");
New_Line;
Result := Sqrt (abs (X));
return Result;
else
Defloat (X, N, F);
Y := Sqrt_C1 + Mantissa_Type (Sqrt_C2 * F);
for J in 1 .. Sqrt_L1 loop
Y := Y / Mantissa_Divisor_2 + Mantissa_Type
((F / Mantissa_Divisor_2) / Y);
end loop;
if (N mod 2) /= 0 then
Y := Mantissa_Type (Sqrt_C3 * Y);
N := N + 1;
end if;
M := N / 2;
Refloat (M, Y, Result);
return Result;
end if;
exception
when others =>
New_Line;
Put (" EXCEPTION IN SQRT, X = ");
raise Program_Error;
--Put (X);
Put (" RETURNED 1.0");
New_Line;
return One;
end Sqrt;
function Cbrt (X : Float) return Float is
M, N : Exponent_Type;
F, Y : Mantissa_Type;
Result : Float;
subtype Index is Integer range 0 ..
100; -- #########################
Cbrt_L1 : Index := 3;
Cbrt_C1 : Mantissa_Type := 0.5874009;
Cbrt_C2 : Mantissa_Type := 0.4125990;
Cbrt_C3 : Mantissa_Type := 0.62996_05249;
Cbrt_C4 : Mantissa_Type := 0.79370_05260;
begin
if X = Zero then
Result := Zero;
return Result;
else
Defloat (X, N, F);
F := abs (F);
Y := Cbrt_C1 + Mantissa_Type (Cbrt_C2 * F);
for J in 1 .. Cbrt_L1 loop
Y := Y - (Y / Mantissa_Divisor_3 -
Mantissa_Type ((F / Mantissa_Divisor_3) /
Mantissa_Type (Y * Y)));
end loop;
case (N mod 3) is
when 0 =>
null;
when 1 =>
Y := Mantissa_Type (Cbrt_C3 * Y);
N := N + 2;
when 2 =>
Y := Mantissa_Type (Cbrt_C4 * Y);
N := N + 1;
when others =>
null;
end case;
M := N / 3;
if X < Zero then
Y := -Y;
end if;
Refloat (M, Y, Result);
return Result;
end if;
exception
when others =>
Result := One;
if X < Zero then
Result := -One;
end if;
New_Line;
Put ("EXCEPTION IN CBRT, X = ");
raise Program_Error;
--Put (X);
Put (" RETURNED ");
raise Program_Error;
--Put (Result);
New_Line;
return Result;
end Cbrt;
function Log (X : Float) return Float is
-- Uses fixed formulation for generality
Result : Float;
N : Exponent_Type;
Xn : Float;
Y : Float;
F : Mantissa_Type;
Z, Zden, Znum : Mantissa_Type;
C0 : constant Mantissa_Type := 0.20710_67811_86547_52440;
-- SQRT(0.5) - 0.5
C1 : constant Float := 8#0.543#;
C2 : constant Float := -2.12194_44005_46905_82767_9E-4;
function R (Z : Mantissa_Type) return Mantissa_Type is
-- Use fixed formulation here because the float coeficents are > 1.0
-- and would exceed the limits on a MANTISSA_TYPE
A0 : constant Mantissa_Type := 0.04862_85276_587;
B0 : constant Mantissa_Type := 0.69735_92187_803;
B1 : constant Mantissa_Type := -0.125;
C : constant Mantissa_Type := 0.01360_09546_862;
begin
return Z + Mantissa_Type
(Z * Mantissa_Type
(Mantissa_Type (Z * Z) *
(C + Mantissa_Type
(A0 / (B0 +
Mantissa_Type
(B1 * Mantissa_Type
(Z * Z)))))));
end R;
begin
if X < Zero then
New_Line;
Put ("CALLED LOG FOR NEGATIVE ");
raise Program_Error; --put (X);;
Put (" USE ABS => ");
Result := Log (abs (X));
raise Program_Error;
--Put (Result);
New_Line;
elsif X = Zero then
New_Line;
Put ("CALLED LOG FOR ZERO ARGUMENT, RETURNED ");
Result := -Xmax; -- SUPPOSED TO BE -LARGE
raise Program_Error;--Put (Result);
New_Line;
else
Defloat (X, N, F);
Znum := F - Mantissa_Half;
Y := Convert_To_Float (Znum);
Zden := Znum / Mantissa_Divisor_2 + Mantissa_Half;
if Znum > C0 then
Y := Y - Mantissa_Half;
Znum := Znum - Mantissa_Half;
Zden := Zden + Mantissa_Half / Mantissa_Divisor_2;
else
N := N - 1;
end if;
Z := Mantissa_Type (Znum / Zden);
Result := Convert_To_Float (R (Z));
if N /= 0 then
Xn := Convert_To_Float (N);
Result := (Xn * C2 + Result) + Xn * C1;
end if;
end if;
return Result;
exception
when others =>
New_Line;
Put (" EXCEPTION IN LOG, X = ");
raise Program_Error; --put (X);;
Put (" RETURNED 0.0");
New_Line;
return Zero;
end Log;
function Log10 (X : Float) return Float is
Log_10_Of_2 : constant Float :=
Convert_To_Float (Mantissa_Type (
8#0.33626_75425_11562_41615#));
begin
return Log (X) * Log_10_Of_2;
end Log10;
function Exp (X : Float) return Float is
Result : Float;
N : Exponent_Type;
Xg, Xn, X1, X2 : Float;
F, G : Mantissa_Type;
Bigx : Float := Exp_Large;
Smallx : Float := Exp_Small;
One_Over_Log_2 : constant Float := 1.4426_95040_88896_34074;
C1 : constant Float := 0.69335_9375;
C2 : constant Float := -2.1219_44400_54690_58277E-4;
function R (G : Mantissa_Type) return Mantissa_Type is
Z, Gp, Q : Mantissa_Type;
P0 : constant Mantissa_Type := 0.24999_99999_9992;
P1 : constant Mantissa_Type := 0.00595_04254_9776;
Q0 : constant Mantissa_Type := 0.5;
Q1 : constant Mantissa_Type := 0.05356_75176_4522;
Q2 : constant Mantissa_Type := 0.00029_72936_3682;
begin
Z := Mantissa_Type (G * G);
Gp := Mantissa_Type ((Mantissa_Type (P1 * Z) + P0) * G);
Q := Mantissa_Type ((Mantissa_Type (Q2 * Z) + Q1) * Z) + Q0;
return Mantissa_Half + Mantissa_Type (Gp / (Q - Gp));
end R;
begin
if X > Bigx then
New_Line;
Put (" EXP CALLED WITH TOO BIG A POSITIVE ARGUMENT, ");
raise Program_Error; --put (X);;
Put (" RETURNED XMAX");
New_Line;
Result := Xmax;
elsif X < Smallx then
New_Line;
Put (" EXP CALLED WITH TOO BIG A NEGATIVE ARGUMENT, ");
raise Program_Error; --put (X);;
Put (" RETURNED ZERO");
New_Line;
Result := Zero;
elsif abs (X) < Eps then
Result := One;
else
N := Exponent_Type (X * One_Over_Log_2);
Xn := Convert_To_Float (N);
X1 := Round (X);
X2 := X - X1;
Xg := ((X1 - Xn * C1) + X2) - Xn * C2;
G := Mantissa_Type (Xg);
N := N + 1;
F := R (G);
Refloat (N, F, Result);
end if;
return Result;
exception
when others =>
New_Line;
Put (" EXCEPTION IN EXP, X = ");
raise Program_Error; --put (X);;
Put (" RETURNED 1.0");
New_Line;
return One;
end Exp;
function "**" (X, Y : Float) return Float is
-- This is the last function to be coded since it appeared that it really
-- was un-Ada-like and ought not be in the regular package
-- Nevertheless it was included in this version
-- It is specific for FLOAT and does not have the MANTISSA_TYPE generality
M, N : Exponent_Type;
G : Mantissa_Type;
P, Temp, Iw1, I : Integer;
Result, Z, V, R, U1, U2, W, W1, W2, W3, Y1, Y2 : Float;
K : constant Float := 0.44269_50408_88963_40736;
Ibigx : constant Integer := Integer
(Truncate (16.0 * Log (Xmax) - 1.0));
Ismallx : constant Integer := Integer
(Truncate (16.0 * Log (Xmin) + 1.0));
P1 : constant Float := 0.83333_32862_45E-1;
P2 : constant Float := 0.12506_48500_52E-1;
Q1 : constant Float := 0.69314_71805_56341;
Q2 : constant Float := 0.24022_65061_44710;
Q3 : constant Float := 0.55504_04881_30765E-1;
Q4 : constant Float := 0.96162_06595_83789E-2;
Q5 : constant Float := 0.13052_55159_42810E-2;
A1 : array (1 .. 17) of Float :=
(8#1.00000_0000#, 8#0.75222_5750#, 8#0.72540_3067#, 8#0.70146_3367#,
8#0.65642_3746#, 8#0.63422_2140#, 8#0.61263_4520#, 8#0.57204_2434#,
8#0.55202_3631#, 8#0.53254_0767#, 8#0.51377_3265#,
8#0.47572_4623#, 8#0.46033_7602#, 8#0.44341_7233#,
8#0.42712_7017#, 8#0.41325_3033#, 8#0.40000_0000#);
A2 : array (1 .. 8) of Float :=
(8#0.00000_00005_22220_66302_61734_72062#,
8#0.00000_00003_02522_47021_04062_61124#,
8#0.00000_00005_21760_44016_17421_53016#,
8#0.00000_00007_65401_41553_72504_02177#,
8#0.00000_00002_44124_12254_31114_01243#,
8#0.00000_00000_11064_10432_66404_42174#,
8#0.00000_00004_72542_16063_30176_55544#,
8#0.00000_00001_74611_03661_23056_22556#);
function Reduce (V : Float) return Float is
begin
return Float (Integer (16.0 * V)) * 0.0625;
end Reduce;
begin
if X <= Zero then
if X < Zero then
Result := (abs (X)) ** Y;
New_Line;
Put ("X**Y CALLED WITH X = ");
raise Program_Error; --put (X);
New_Line;
Put ("USED ABS, RETURNED ");
raise Program_Error;--Put (Result);
New_Line;
else
if Y <= Zero then
if Y = Zero then
Result := Zero;
else
Result := Xmax;
end if;
New_Line;
Put ("X**Y CALLED WITH X = 0, Y = ");
raise Program_Error; --put (Y);
New_Line;
Put ("RETURNED ");
raise Program_Error;--Put (Result);
New_Line;
else
Result := Zero;
end if;
end if;
else
Defloat (X, M, G);
P := 1;
if G <= A1 (9) then
P := 9;
end if;
if G <= A1 (P + 4) then
P := P + 4;
end if;
if G <= A1 (P + 2) then
P := P + 2;
end if;
Z := ((G - A1 (P + 1)) - A2 ((P + 1) / 2)) / (G + A1 (P + 1));
Z := Z + Z;
V := Z * Z;
R := (P2 * V + P1) * V * Z;
R := R + K * R;
U2 := (R + Z * K) + Z;
U1 := Float (Integer (M) * 16 - P) * 0.0625;
Y1 := Reduce (Y);
Y2 := Y - Y1;
W := U2 * Y + U1 * Y2;
W1 := Reduce (W);
W2 := W - W1;
W := W1 + U1 * Y1;
W1 := Reduce (W);
W2 := W2 + (W - W1);
W3 := Reduce (W2);
Iw1 := Integer (Truncate (16.0 * (W1 + W3)));
W2 := W2 - W3;
if W > Float (Ibigx) then
Result := Xmax;
Put ("X**Y CALLED X =");
raise Program_Error; --put (X);
Put (" Y =");
raise Program_Error; --put (Y);
Put (" TOO LARGE RETURNED ");
raise Program_Error;--Put (Result);
New_Line;
elsif W < Float (Ismallx) then
Result := Zero;
Put ("X**Y CALLED X =");
raise Program_Error;
--put(X);
Put (" Y =");
raise Program_Error;
Put (" TOO SMALL RETURNED ");
raise Program_Error;--Put (Result);
New_Line;
else
if W2 > Zero then
W2 := W2 - 0.0625;
Iw1 := Iw1 + 1;
end if;
if Iw1 < Integer (Zero) then
I := 0;
else
I := 1;
end if;
M := Exponent_Type (I + Iw1 / 16);
P := 16 * Integer (M) - Iw1;
Z := ((((Q5 * W2 + Q4) * W2 + Q3) * W2 + Q2) * W2 + Q1) * W2;
Z := A1 (P + 1) + (A1 (P + 1) * Z);
Refloat (M, Z, Result);
end if;
end if;
return Result;
end "**";
begin
Exp_Large := Log (Xmax) * (One - Eps);
Exp_Small := Log (Xmin) * (One - Eps);
end Core_Functions;