|
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: 18727 (0x4927)
Types: TextFile
Names: »B«
└─⟦180fe333a⟧ Bits:30000405 8mm tape, Rational 1000, SW CATALOG, 10_20_0
└─⟦180fe333a⟧ Bits:30000537 8mm tape, Rational 1000, SW Catalog 10_20_0
└─⟦5cb1d1d7f⟧ »DATA«
└─⟦3b1ee7bd8⟧
└─⟦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;
with Core_Functions;
use Core_Functions;
package body Trig_Lib is
-- PRELIMINARY VERSION *********************************
-- 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
-- This particular version is stripped to work with FLOAT and INTEGER
-- and uses a mantissa represented as a FLOAT
-- 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
-- 16 JULY 1982 W A WHITAKER AFATL EGLIN AFB FL 32542
-- T C EICHOLTZ USAFA
function Sin (X : Float) return Float is
Sgn, Y : Float;
N : Integer;
Xn : Float;
F, G, X1, X2 : Float;
Result : Float;
Ymax : Float := Float (Integer (Pi * Two ** (It / 2)));
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
C1 : constant Float := 3.140625;
C2 : constant Float := 9.6765_35897_93E-4;
function R (G : Float) return Float is
R1 : constant Float := -0.16666_66660_883;
R2 : constant Float := 0.83333_30720_556E-2;
R3 : constant Float := -0.19840_83282_313E-3;
R4 : constant Float := 0.27523_97106_775E-5;
R5 : constant Float := -0.23868_34640_601E-7;
begin
return ((((R5 * G + R4) * G + R3) * G + R2) * G + R1) * G;
end R;
begin
if X < Zero then
Sgn := -One;
Y := -X;
else
Sgn := One;
Y := X;
end if;
if Y > Ymax then
New_Line;
Put (" SIN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
Put (X);
New_Line;
end if;
N := Integer (Y * One_Over_Pi);
Xn := Convert_To_Float (N);
if N mod 2 /= 0 then
Sgn := -Sgn;
end if;
X1 := Truncate (abs (X));
X2 := abs (X) - X1;
F := ((X1 - Xn * C1) + X2) - Xn * C2;
if abs (F) < Epsilon then
Result := F;
else
G := F * F;
Result := F + F * R (G);
end if;
return (Sgn * Result);
end Sin;
function Cos (X : Float) return Float is
Sgn, Y : Float;
N : Integer;
Xn : Float;
F, G, X1, X2 : Float;
Result : Float;
Ymax : Float := Float (Integer (Pi * Two ** (It / 2)));
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
C1 : constant Float := 3.140625;
C2 : constant Float := 9.6765_35897_93E-4;
function R (G : Float) return Float is
R1 : constant Float := -0.16666_66660_883;
R2 : constant Float := 0.83333_30720_556E-2;
R3 : constant Float := -0.19840_83282_313E-3;
R4 : constant Float := 0.27523_97106_775E-5;
R5 : constant Float := -0.23868_34640_601E-7;
begin
return ((((R5 * G + R4) * G + R3) * G + R2) * G + R1) * G;
end R;
begin
Sgn := 1.0;
Y := abs (X) + Pi_Over_Two;
if Y > Ymax then
New_Line;
Put (" COS CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
Put (X);
New_Line;
end if;
N := Integer (Y * One_Over_Pi);
Xn := Convert_To_Float (N);
if N mod 2 /= 0 then
Sgn := -Sgn;
end if;
Xn := Xn - 0.5; -- TO FORM COS INSTEAD OF SIN
X1 := Truncate (abs (X));
X2 := abs (X) - X1;
F := ((X1 - Xn * C1) + X2) - Xn * C2;
if abs (F) < Epsilon then
Result := F;
else
G := F * F;
Result := F + F * R (G);
end if;
return (Sgn * Result);
end Cos;
function Tan (X : Float) return Float is
Sgn, Y : Float;
N : Integer;
Xn : Float;
F, G, X1, X2 : Float;
Result : Float;
Ymax : Float := Float (Integer (Pi * Two ** (It / 2))) / 2.0;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
C1 : constant Float := 8#1.444#;
C2 : constant Float := 4.8382_67948_97E-4;
function R (G : Float) return Float is
P0 : constant Float := 1.0;
P1 : constant Float := -0.11136_14403_566;
P2 : constant Float := 0.10751_54738_488E-2;
Q0 : constant Float := 1.0;
Q1 : constant Float := -0.44469_47720_281;
Q2 : constant Float := 0.15973_39213_300E-1;
begin
return ((P2 * G + P1) * G * F + F) /
(((Q2 * G + Q1) * G + 0.5) + 0.5);
end R;
begin
Y := abs (X);
if Y > Ymax then
New_Line;
Put (" TAN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
Put (X);
New_Line;
end if;
N := Integer (X * Two_Over_Pi);
Xn := Convert_To_Float (N);
X1 := Truncate (X);
X2 := X - X1;
F := ((X1 - Xn * C1) + X2) - Xn * C2;
if abs (F) < Epsilon then
Result := F;
else
G := F * F;
Result := R (G);
end if;
if N mod 2 = 0 then
return Result;
else
return -1.0 / Result;
end if;
end Tan;
function Cot (X : Float) return Float is
Sgn, Y : Float;
N : Integer;
Xn : Float;
F, G, X1, X2 : Float;
Result : Float;
Ymax : Float := Float (Integer (Pi * Two ** (It / 2))) / 2.0;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
Epsilon1 : Float := 1.0 / Xmax;
C1 : constant Float := 8#1.444#;
C2 : constant Float := 4.8382_67948_97E-4;
function R (G : Float) return Float is
P0 : constant Float := 1.0;
P1 : constant Float := -0.11136_14403_566;
P2 : constant Float := 0.10751_54738_488E-2;
Q0 : constant Float := 1.0;
Q1 : constant Float := -0.44469_47720_281;
Q2 : constant Float := 0.15973_39213_300E-1;
begin
return ((P2 * G + P1) * G * F + F) /
(((Q2 * G + Q1) * G + 0.5) + 0.5);
end R;
begin
Y := abs (X);
if Y < Epsilon1 then
New_Line;
Put (" COT CALLED WITH ARGUMENT TOO NEAR ZERO ");
Put (X);
New_Line;
if X < 0.0 then
return -Xmax;
else
return Xmax;
end if;
end if;
if Y > Ymax then
New_Line;
Put (" COT CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
Put (X);
New_Line;
end if;
N := Integer (X * Two_Over_Pi);
Xn := Convert_To_Float (N);
X1 := Truncate (X);
X2 := X - X1;
F := ((X1 - Xn * C1) + X2) - Xn * C2;
if abs (F) < Epsilon then
Result := F;
else
G := F * F;
Result := R (G);
end if;
if N mod 2 /= 0 then
return -Result;
else
return 1.0 / Result;
end if;
end Cot;
function Asin (X : Float) return Float is
G, Y : Float;
Result : Float;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
function R (G : Float) return Float is
P1 : constant Float := -0.27516_55529_0596E1;
P2 : constant Float := 0.29058_76237_4859E1;
P3 : constant Float := -0.59450_14419_3246;
Q0 : constant Float := -0.16509_93320_2424E2;
Q1 : constant Float := 0.24864_72896_9164E2;
Q2 : constant Float := -0.10333_86707_2113E2;
Q3 : constant Float := 1.0;
begin
return (((P3 * G + P2) * G + P1) * G) /
(((G + Q2) * G + Q1) * G + Q0);
end R;
begin
Y := abs (X);
if Y > Half then
if Y > 1.0 then
New_Line;
Put (" ASIN CALLED FOR ");
Put (X);
Put (" (> 1) TRUNCATED TO 1, CONTINUED");
New_Line;
Y := 1.0;
end if;
G := ((0.5 - Y) + 0.5) / 2.0;
Y := -2.0 * Sqrt (G);
Result := Y + Y * R (G);
Result := (Pi_Over_Four + Result) + Pi_Over_Four;
else
if Y < Epsilon then
Result := Y;
else
G := Y * Y;
Result := Y + Y * R (G);
end if;
end if;
if X < 0.0 then
Result := -Result;
end if;
return Result;
end Asin;
function Acos (X : Float) return Float is
G, Y : Float;
Result : Float;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
function R (G : Float) return Float is
P1 : constant Float := -0.27516_55529_0596E1;
P2 : constant Float := 0.29058_76237_4859E1;
P3 : constant Float := -0.59450_14419_3246;
Q0 : constant Float := -0.16509_93320_2424E2;
Q1 : constant Float := 0.24864_72896_9164E2;
Q2 : constant Float := -0.10333_86707_2113E2;
Q3 : constant Float := 1.0;
begin
return (((P3 * G + P2) * G + P1) * G) /
(((G + Q2) * G + Q1) * G + Q0);
end R;
begin
Y := abs (X);
if Y > Half then
if Y > 1.0 then
New_Line;
Put (" ACOS CALLED FOR ");
Put (X);
Put (" (> 1) TRUNCATED TO 1, CONTINUED");
New_Line;
Y := 1.0;
end if;
G := ((0.5 - Y) + 0.5) / 2.0;
Y := -2.0 * Sqrt (G);
Result := Y + Y * R (G);
if X < 0.0 then
Result := (Pi_Over_Two + Result) + Pi_Over_Two;
else
Result := -Result;
end if;
else
if Y < Epsilon then
Result := Y;
else
G := Y * Y;
Result := Y + Y * R (G);
end if;
if X < 0.0 then
Result := (Pi_Over_Four + Result) + Pi_Over_Four;
else
Result := (Pi_Over_Four - Result) + Pi_Over_Four;
end if;
end if;
return Result;
end Acos;
function Atan (X : Float) return Float is
F, G : Float;
subtype Region is Integer range 0 .. 3; -- ##########
N : Region;
Result : Float;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
Sqrt_3 : constant Float := 1.73205_08075_68877_29353;
Sqrt_3_Minus_1 : constant Float := 0.73205_08075_68877_29353;
Two_Minus_Sqrt_3 : constant Float := 0.26794_91924_31122_70647;
function R (G : Float) return Float is
P0 : constant Float := -0.14400_83448_74E1;
P1 : constant Float := -0.72002_68488_98;
Q0 : constant Float := 0.43202_50389_19E1;
Q1 : constant Float := 0.47522_25845_99E1;
Q2 : constant Float := 1.0;
begin
return ((P1 * G + P0) * G) / ((G + Q1) * G + Q0);
end R;
begin
F := abs (X);
if F > 1.0 then
F := 1.0 / F;
N := 2;
else
N := 0;
end if;
if F > Two_Minus_Sqrt_3 then
F := (((Sqrt_3_Minus_1 * F - 0.5) - 0.5) + F) / (Sqrt_3 + F);
N := N + 1;
end if;
if abs (F) < Epsilon then
Result := F;
else
G := F * F;
Result := F + F * R (G);
end if;
if N > 1 then
Result := -Result;
end if;
case N is
when 0 =>
Result := Result;
when 1 =>
Result := Pi_Over_Six + Result;
when 2 =>
Result := Pi_Over_Two + Result;
when 3 =>
Result := Pi_Over_Three + Result;
end case;
if X < 0.0 then
Result := -Result;
end if;
return Result;
end Atan;
function Atan2 (V, U : Float) return Float is
X, Result : Float;
begin
if U = 0.0 then
if V = 0.0 then
Result := 0.0;
New_Line;
Put (" ATAN2 CALLED WITH 0/0 RETURNED ");
Put (Result);
New_Line;
elsif V > 0.0 then
Result := Pi_Over_Two;
else
Result := -Pi_Over_Two;
end if;
else
X := abs (V / U);
-- If underflow or overflow is detected, go to the exception
Result := Atan (X);
if U < 0.0 then
Result := Pi - Result;
end if;
if V < 0.0 then
Result := -Result;
end if;
end if;
return Result;
exception
when Numeric_Error =>
if abs (V) > abs (U) then
Result := Pi_Over_Two;
if V < 0.0 then
Result := -Result;
end if;
else
Result := 0.0;
if U < 0.0 then
Result := Pi - Result;
end if;
end if;
return Result;
end Atan2;
function Sinh (X : Float) return Float is
G, W, Y, Z : Float;
Result : Float;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
Ybar : Float := Exp_Large;
Ln_V : Float := 8#0.542714#;
V_Over_2_Minus_1 : Float := 0.13830_27787_96019_02638E-4;
Wmax : Float := Ybar - Ln_V + 0.69;
function R (G : Float) return Float is
P0 : constant Float := 0.10622_28883_7151E4;
P1 : constant Float := 0.31359_75645_6058E2;
P2 : constant Float := 0.34364_14035_8506;
Q0 : constant Float := 0.63733_73302_1822E4;
Q1 : constant Float := -0.13051_01250_9199E3;
Q2 : constant Float := 1.0;
begin
return (((P2 * G + P1) * G + P0) * G) / ((G + Q1) * G + Q0);
end R;
begin
Y := abs (X);
if Y <= 1.0 then
if Y < Epsilon then
Result := X;
else
G := X * X;
Result := X + X * R (G);
end if;
else
if Y <= Ybar then
Z := Exp (Y);
Result := (Z - 1.0 / Z) / 2.0;
else
W := Y - Ln_V;
if W > Wmax then
New_Line;
Put (" SINH CALLED WITH TOO LARGE ARGUMENT ");
Put (X);
Put (" RETURN BIG");
New_Line;
W := Wmax;
end if;
Z := Exp (W);
Result := Z + V_Over_2_Minus_1 * Z;
end if;
if X < 0.0 then
Result := -Result;
end if;
end if;
return Result;
end Sinh;
function Cosh (X : Float) return Float is
G, W, Y, Z : Float;
Result : Float;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
Ybar : Float := Exp_Large;
Ln_V : Float := 8#0.542714#;
V_Over_2_Minus_1 : Float := 0.13830_27787_96019_02638E-4;
Wmax : Float := Ybar - Ln_V + 0.69;
function R (G : Float) return Float is
P0 : constant Float := 0.10622_28883_7151E4;
P1 : constant Float := 0.31359_75645_6058E2;
P2 : constant Float := 0.34364_14035_8506;
Q0 : constant Float := 0.63733_73302_1822E4;
Q1 : constant Float := -0.13051_01250_9199E3;
Q2 : constant Float := 1.0;
begin
return (((P2 * G + P1) * G + P0) * G) / ((G + Q1) * G + Q0);
end R;
begin
Y := abs (X);
if Y <= Ybar then
Z := Exp (Y);
Result := (Z + 1.0 / Z) / 2.0;
else
W := Y - Ln_V;
if W > Wmax then
New_Line;
Put (" COSH CALLED WITH TOO LARGE ARGUMENT ");
Put (X);
Put (" RETURN BIG");
New_Line;
W := Wmax;
end if;
Z := Exp (W);
Result := Z + V_Over_2_Minus_1 * Z;
end if;
return Result;
end Cosh;
function Tanh (X : Float) return Float is
G, W, Y, Z : Float;
Result : Float;
Beta : Float := Convert_To_Float (Ibeta);
Epsilon : Float := Beta ** (-It / 2);
Xbig : Float :=
(Log (2.0) + Convert_To_Float (It + 1) * Log (Beta)) / 2.0;
Ln_3_Over_2 : Float := 0.54930_61443_34054_84570;
function R (G : Float) return Float is
P0 : constant Float := -0.21063_95800_0245E2;
P1 : constant Float := -0.93363_47565_2401;
Q0 : constant Float := 0.63191_87401_5582E2;
Q1 : constant Float := 0.28077_65347_0471E2;
Q2 : constant Float := 1.0;
begin
return ((P1 * G + P0) * G) / ((G + Q1) * G + Q0);
end R;
begin
Y := abs (X);
if Y > Xbig then
Result := 1.0;
else
if Y > Ln_3_Over_2 then
Result := 0.5 - 1.0 / (Exp (Y + Y) + 1.0);
Result := Result + Result;
else
if Y < Epsilon then
Result := Y;
else
G := Y * Y;
Result := Y + Y * R (G);
end if;
end if;
end if;
if X < 0.0 then
Result := -Result;
end if;
return Result;
end Tanh;
begin
null;
end Trig_Lib;