|
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: 57344 (0xe000)
Types: Ada Source
Notes: 03_class, FILE, R1k_Segment, e3_tag, package body Generic_Elementary_Functions, seg_0069cf
└─⟦8527c1e9b⟧ Bits:30000544 8mm tape, Rational 1000, Arrival backup of disks in PAM's R1000
└─⟦5a81ac88f⟧ »Space Info Vol 1«
└─⟦this⟧
-- This package body contains an implementation of the proposed standard for
-- elementary mathematical functions, as defined by the ISO-IEC/JTC1/SC22/WG9
-- (Ada) committee (draft 1.2, December 1990)
--
-- Users recognize that the Math Support packages provided as part of the
-- Rational Environment are delivered "AS IS". No warranties are given,
-- whether express, implied, or statutory, including implied warranties of
-- fitness for a particular purpose or mechantability. In no event will
-- Rational be liable in tort, negligence, or other liability incurred as a
-- result of the use of this Math Support packages.
--
-- In addition, Rational does not guarantee conformance to the proposed
-- standards, especially in the accuracy of the results of the various math
-- functions.
--
-- This code can be distributed and modified freely, although Rational
-- would be happy to receive comments, bug reports, and suggestions for
-- improvements at:
-- Rational
-- 3320 Scott Boulevard
-- Santa Clara, CA 95054 USA
-- or by e-mail at:
-- Support@Rational.COM
--
with Primitive_Functions;
package body Generic_Elementary_Functions is
Pi : constant :=
3.1415926535897932384626433832795028841971693993750990042;
Pi_Over_2 : constant := Pi / 2.0;
Pi_Over_4 : constant := Pi / 4.0;
Two_Pi : constant := Pi * 2.0;
One_Over_Pi : constant := 1.0 / Pi;
One_Over_Two_Pi : constant := 1.0 / Two_Pi;
Two_Over_Pi : constant := 2.0 / Pi;
High_Threshold : constant Float :=
Float (Float'Machine_Radix) ** (Float'Machine_Mantissa / 2);
Low_Threshold : constant Float :=
Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa / 2);
Log_Float_Large : constant := 141.402024;-- Log(float'large);
Log_2 : constant Float :=
0.693147180559945309417232121458176568075500133107812;
-- constants for "**"
--
type Float_Array is array (Natural range <>) of Float;
Expx_A1 : Float_Array (1 .. 17);
Expx_A2 : Float_Array (1 .. 8);
Sqrt_One_Half : constant Float :=
0.7071067811865475244008443621048490392848359375608;
------------------------------------------------------------------------------
use Primitive_Functions;
function Is_Odd (F : Float) return Boolean is
-- assume F is already an integer !!
F_Over_2 : Float := F / 2.0;
begin
return (Truncate (F_Over_2) /= F_Over_2);
end Is_Odd;
-------------------------------------------------------------------------------
function Sqrt (X : Float_Type) return Float_Type is
-- Cody & Waite, pp.17-34
-- Hart & al., pp.89-96
A : Float := Float (X);
Adj : Float;
F : Float;
N : Integer;
Y : Float;
Bits : Natural;
begin
if X < 0.0 then
raise Argument_Error;
elsif X = 0.0 then
return 0.0;
end if;
----Break A into a fraction and an exponent.
Decompose (A, F, N);
----Compute an approximation of the answer for sqrt(f).
Y := 0.41731 + 0.59016 * F;
----Now we put back the exponential information.
if N rem 2 /= 0 then
Y := Y * 8#0.55202_36314_77473_63110#;
N := N + 1;
end if;
N := N / 2;
Y := Y * Float (Float'Machine_Radix) ** N;
----Now use Newton's approximation to compute the result of sqrt(A).
Bits := 8; -- Currently valid bits (minimum)
while Bits < Float'Size loop
Y := 0.5 * (Y + A / Y);
Bits := Bits + Bits;
end loop;
----We now have an answer that is correct to within the bottom bit (or two).
-- See if we can improve on that. The Newton's approximation cannot
-- do any better because the hardware does (or does not do) rounding and this
-- introduces a "granularity" to the continuous Newton approximation.
-- So, we will sit and diddle with the bottom bit of the answer. We want a
-- value that when squared is strictly less than our original argument.
F := Y * Y;
if F /= A then
Adj := Float (Float'Machine_Radix) ** (N - Float'Machine_Mantissa);
----Adj is so small that it only diddles the bottom bit of the
-- mantissa; we so assert.
begin
while F < A loop
Y := Y + Adj;
F := Y * Y;
end loop;
exception
when Numeric_Error | Constraint_Error =>
----We got an agument that was at/close-to Float'Large
-- Reduce Y back by Adj and return it.
Y := Y - Adj;
return Float_Type (Y);
end;
while F > A loop
Y := Y - Adj;
F := Y * Y;
end loop;
end if;
----Return the correct result.
return Float_Type (Y);
end Sqrt;
function Log (X : Float_Type) return Float_Type is
-- Cody & Waite, pp.35-48
A : Float := Float (X);
F : Float;
N : Integer;
Znum : Float;
Zden : Float;
begin
if X < 0.0 then
raise Argument_Error;
end if;
if X = 0.0 then
raise Constraint_Error;
end if;
----Extract the fraction and the exponent from the argument. Adjust them
-- so that 0.5 <= f < 1.0. We want the exponent expressed as a power of
-- 2.0.
Decompose (A, F, N);
----Now we adjust N and z if f <= sqrt(0.5) and we just adjust z if
-- f > sqrt(0.5).
if F > Sqrt_One_Half then
Znum := (F - 0.5) - 0.5;
Zden := F * 0.5 + 0.5;
else
N := N - 1;
Znum := F - 0.5;
Zden := Znum * 0.5 + 0.5;
end if;
----Now we evaluate a rational approximation of exp(z).
Znum := Znum / Zden;
F := Znum * Znum;
Zden :=
((-0.78956_11288_74912_57267E+00 * F + 0.16383_94356_30215_34222E+02)
* F - 0.64124_94342_37455_81147E+02)
/ (((F - 0.35667_97773_90346_46171E+02)
* F + 0.31203_22209_19245_32844E+03)
* F - 0.76949_93210_84948_79777E+03);
F := Znum + Znum * (F * Zden);
----Now compute the result.
Znum := Float (N);
F := (-2.12194_44005_46905_827679E-04 * Znum + F);
F := F + Znum * 8#0.543#;
return Float_Type (F);
end Log;
function Log (X, Base : Float_Type) return Float_Type is
begin
if X < 0.0 or Base <= 0.0 or Base = 1.0 then
raise Argument_Error;
end if;
if X = 0.0 then
raise Constraint_Error;
end if;
return Log (X) / Log (Base);
end Log;
function Exp (X : Float_Type) return Float_Type is
-- Hart et. al., pp.97-105
A : Float := Float (X);
N : Integer;
Xn : Float;
G : Float;
Xi : Integer;
X1 : Float;
X2 : Float;
Pz : Float;
Qz : Float;
Ln_Float_Large : constant := 141.402024;-- Log(float'large);
Ln_Float_Small : constant := -711.169; -- log(float'safe_small);
begin
----Make sure that it is possible to represent the result. (too big/small)
if A > Ln_Float_Large then
raise Constraint_Error;
end if;
if A < Ln_Float_Small then
return 0.0;
end if;
----If the answer will inevitably be 1.0 then return 1.0.
if abs A <
Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa) / 2.0 then
return 1.0;
end if;
----Perform argument range reduction.
N := Integer (A * 1.44269504088896340735992468100189213742664595602385);
Xn := Float (N);
Xi := Integer (abs A);
if Float (Xi) > abs A then
Xi := Xi - 1;
end if;
if A < 0.0 then
X1 := -Float (Xi);
else
X1 := Float (Xi);
end if;
X2 := A - X1;
G := ((X1 - Xn * 8#0.543#) + X2) + Xn * 2.1219_44400_54690_58277E-04;
----Now approximate exp(g)
X1 := G * G;
Pz := (0.16520_33002_68279_130E-04 * X1 +
0.69436_00015_11792_852E-02) * X1 + 0.24999_99999_99999_993E+00;
Qz := (0.49586_28849_05441_294E-03 * X1 +
0.55553_86669_69001_118E-01) * X1 + 0.5;
Pz := Pz * G;
G := 0.5 + Pz / (Qz - Pz);
----Compute the final result.
N := N + 1;
Xi := N / 2;
N := N - Xi;
G := G * 2.0 ** Xi;
G := G * 2.0 ** N;
return Float_Type (G);
end Exp;
function "**" (X, Y : Float_Type) return Float_Type is
M : Integer;
F : Float;
P : Integer;
Z : Float;
Rv : Float;
U1 : Float;
U2 : Float; Y1 : Float;
Y2 : Float;
Iw1 : Integer;
W : Float;
W1 : Float;
W2 : Float;
I : Integer;
begin
if X = 0.0 and Y /= 0.0 then
raise Constraint_Error;
end if;
if X < 0.0 then
raise Argument_Error;
end if;
if Y = 0.0 then
return 1.0;
end if;
Decompose (Float (X), F, M);
----Compute the P index.
P := 1;
if F <= Expx_A1 (9) then
P := 9;
end if;
if F <= Expx_A1 (P + 4) then
P := P + 4;
end if;
if F <= Expx_A1 (P + 2) then
P := P + 2;
end if;
----Compute z.
Z := F - Expx_A1 (P + 1);
Z := Z - Expx_A2 ((P + 1) / 2);
Z := Z / (F + Expx_A1 (P + 1));
Z := Z + Z;
----Now compute the z*P(z**2) approximation.
F := Z * Z;
Rv :=
((0.43445_77567_21631_19635E-03 * F + 0.22321_42128_59242_58967E-02) *
F + 0.12500_00000_05037_99174E-01) * F +
0.83333_33333_33332_11405E-01;
Rv := Z * F * Rv;
Rv := Rv + 0.44269_50408_88963_40736 * Rv;
U2 := (Rv + Z * 0.44269_50408_88963_40736) + Z;
----Now create W, our partial answer.
U1 := Float (M * 16 - P) * 0.0625;
Y1 := Truncate (Float (Y) * 16.0) * 0.0625;
Y2 := Float (Y) - Y1;
W := U2 * Float (Y) + U1 * Y2;
W1 := Truncate (W * 16.0) * 0.0625;
W2 := W - W1;
W := W1 + U1 * Y1;
W1 := Truncate (W * 16.0) * 0.0625;
W2 := W2 + (W - W1);
W := Truncate (W2 * 16.0) * 0.0625;
Iw1 := Integer (Truncate (16.0 * (W1 + W)));
W2 := W2 - W;
if W2 > 0.0 then
Iw1 := Iw1 + 1;
W2 := W2 - 1.0 / 16.0;
end if;
if Iw1 < 0 then
I := 0;
else
I := 1;
end if;
----Compute the new p,r,m values.
M := Iw1 / 16 + I;
P := 16 * M - Iw1;
----Now do Z.
Z := ((((((0.14928_85268_05956_08186E-04 * W2 +
0.15400_29044_09897_64601E-03) * W2 +
0.13333_54131_35857_84703E-02) * W2 +
0.96181_29059_51724_16964E-02) * W2 +
0.55504_10866_40855_95326E-01) * W2 +
0.24022_65069_59095_37056E+00) * W2 +
0.69314_71805_59945_29629E+00) * W2;
----Nearly there; one more step.
Z := Expx_A1 (P + 1) + Expx_A1 (P + 1) * Z;
I := M / 2;
M := M - I;
Z := Z * Float (Float'Machine_Radix) ** I;
Z := Z * Float (Float'Machine_Radix) ** M;
return Float_Type (Z);
end "**";
function Sin (X : Float_Type) return Float_Type is
R : Float;
F : Float;
G : Float;
Y : Float := Float (X);
Xn : Float;
Tmp : Float;
Sign : Boolean;
begin
if Y < 0.0 then
Y := -Y;
Sign := True;
else
Sign := False;
end if;
----Check for ridiculously large values of Y.
if Y >= Pi * High_Threshold then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (Y * One_Over_Pi);
if Is_Odd (Xn) then
Sign := not Sign;
end if;
----Compute the f portion of the argument.
Tmp := Truncate (Y);
F := ((Tmp - Xn * 8#3.1104#) + (Y - Tmp)) +
Xn * 8.9089_10206_76153_73566_17E-06;
----See if there is no point to computing sin(x).
if abs F < Low_Threshold then
if Sign then
F := -F;
end if;
return Float_Type (F);
end if;
----Compute an answer.
G := F * F;
R := (((((((0.27204_79095_78888_46175E-014 * G -
0.76429_17806_89104_67734E-012) * G +
0.16058_93649_03715_89114E-09) * G -
0.25052_10679_82745_84544E-07) * G +
0.27557_31921_01527_56119E-05) * G -
0.19841_26984_12018_40457E-03) * G +
0.83333_33333_33316_50314E-02) * G -
0.16666_66666_66666_65052E+00) * G;
----Now finish off the result.
R := F + F * R;
if Sign then
R := -R;
end if;
return Float_Type (R);
end Sin;
function Sin (X, Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Sin ((Two_Pi * X) / Cycle);
end Sin;
function Cos (X : Float_Type) return Float_Type is
R : Float;
F : Float;
G : Float;
Y : Float;
Xn : Float;
Tmp : Float;
Sign : Boolean;
begin
Y := Float (abs X) + 1.57079_63267_94896_61923;
Sign := False;
----Check for ridiculously large values of Y.
if Y >= Pi * High_Threshold then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (Y * One_Over_Pi);
if Is_Odd (Xn) then
Sign := not Sign;
end if;
Xn := Xn - 0.5;
----Compute the f portion of the argument.
Tmp := Truncate (Float (abs X));
F := ((Tmp - Xn * 8#3.1104#) + (Float (abs X) - Tmp)) +
Xn * 8.9089_10206_76153_73566_17E-06;
----See if there is no point to computing cos(x).
if abs F < Low_Threshold then
if Sign then
F := -F;
end if;
return Float_Type (F);
end if;
----Compute an answer.
G := F * F;
R := (((((((0.27204_79095_78888_46175E-014 * G -
0.76429_17806_89104_67734E-012) * G +
0.16058_93649_03715_89114E-09) * G -
0.25052_10679_82745_84544E-07) * G +
0.27557_31921_01527_56119E-05) * G -
0.19841_26984_12018_40457E-03) * G +
0.83333_33333_33316_50314E-02) * G -
0.16666_66666_66666_65052E+00) * G;
----Now finish off the result.
R := F + F * R;
if Sign then
R := -R;
end if;
return Float_Type (R);
end Cos;
function Cos (X, Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Cos ((Two_Pi * X) / Cycle);
end Cos;
function Tan (X : Float_Type) return Float_Type is
A : Float := Float (X);
F : Float;
X1 : Float;
Xnum : Float;
Xden : Float;
Even : Boolean;
Xn : Float;
begin
F := abs A;
----See if our argument is much too large to make any sense.
if F > Truncate (High_Threshold * Pi_Over_2) then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (A * Two_Over_Pi);
Even := not Is_Odd (Xn);
----Compute f.
X1 := Truncate (A);
F := ((X1 - Xn * 1.57079) + (A - X1)) +
Xn * 6.32679489661923132169163975144E-06;
----See if f is very very small.
if abs F < Low_Threshold then
Xnum := F;
Xden := 1.0;
else
X1 := F * F;
Xnum := (((-0.17861_70734_22544_26711E-04 * X1 +
0.34248_87823_58905_89960E-02) * X1 -
0.13338_35000_64219_60681E+00) * X1) * F + F;
Xden := (((0.49819_43399_37865_12270E-06 * X1 -
0.31181_53190_70100_27307E-03) * X1 +
0.25663_83228_94401_12864E-01) * X1 -
0.46671_68333_97552_94240E+00) * X1 +
0.5;
Xden := Xden + 0.5;
end if;
----If XN was odd then change sign of result.
if not Even then
Xnum := Xden / (-Xnum);
else
Xnum := Xnum / Xden;
end if;
return Float_Type (Xnum);
end Tan;
function Tan (X, Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Tan (Two_Pi * X / Cycle);
end Tan;
function Cot (X : Float_Type) return Float_Type is
A : Float := Float (X);
F : Float;
X1 : Float;
Xnum : Float;
Xden : Float;
Even : Boolean;
Xn : Float;
----Cot_Eps1 is smallest positive number that can be reciprocated without
-- Numeric_Error.
Cot_Eps1 : constant Float :=
Float (Float'Machine_Radix) ** (Float'Machine_Emin - 1);
begin
F := abs A;
----See if our argument is much too small to handle.
if F < Cot_Eps1 then
if A = 0.0 then
raise Constraint_Error;
elsif A < 0.0 then
return Float_Type (Float'First);
else
return Float_Type (Float'Last);
end if;
end if;
----See if our argument is much too large to make any sense.
if F > Truncate (High_Threshold * Pi_Over_2) then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (A * Two_Over_Pi);
Even := not Is_Odd (Xn);
----Compute f.
X1 := Truncate (A);
F := ((X1 - Xn * 8#1.4442#) + (A - X1)) +
Xn * 4.4544_55103_38706_86783_08E-06;
----See if f is very very small.
if abs F < Low_Threshold then
Xnum := F;
Xden := 1.0;
else
X1 := F * F;
Xnum := (((-0.17861_70734_22544_26711E-04 * X1 +
0.34248_87823_58905_89960E-02) * X1 -
0.13338_35000_64219_60681E+00) * X1) * F + F;
Xden := (((0.49819_43399_37865_12270E-06 * X1 -
0.31181_53190_70100_27307E-03) * X1 +
0.25663_83228_94401_12864E-01) * X1 -
0.46671_68333_97552_94240E+00) * X1 +
0.5;
Xden := Xden + 0.5;
end if;
----If XN was odd then change sign of result.
if not Even then
Xnum := -Xnum / Xden;
else
Xnum := Xden / Xnum;
end if;
return Float_Type (Xnum);
end Cot;
function Cot (X, Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Cot (Two_Pi * X / Cycle);
end Cot;
function Arcsin (X : Float_Type) return Float_Type is
Y : Float := Float (abs X);
F : Float;
G : Float;
I : Integer;
begin
if abs X > 1.0 then
raise Argument_Error;
end if;
----Very small arguments are equal to the answer.
if Y < Low_Threshold then
I := 0;
G := Y;
else
----Compute the intermediate G value differently depending upon the
-- magnitude of the argument.
if Y > 0.5 then
I := 1;
F := ((0.5 - Y) + 0.5) / 2.0;
Y := Float (Sqrt (Float_Type (F)));
Y := -(Y + Y);
G := F;
else
I := 0;
G := Y * Y;
end if;
----Approximate our result.
F := ((((-0.69674_57344_73506_46411E+00 * G +
0.10152_52223_38064_63645E+02) * G -
0.39688_86299_75048_77339E+02) * G +
0.57208_22787_78917_31407E+02) * G -
0.27368_49452_41642_55994E+02) * G
/ (((((G - 0.23823_85915_36702_38830E+02) * G +
0.15095_27084_10306_04719E+03) * G -
0.38186_30336_17501_49284E+03) * G +
0.41714_43024_82604_12556E+03) * G -
0.16421_09671_44985_60795E+03);
G := Y + Y * F;
----Adjust our answer to a different quadrant if necessary.
if I = 1 then
G := 8#0.62207_73250_42055_06043# +
(8#0.62207_73250_42055_06043# + G);
end if;
end if;
----Return the result with the correct sign.
if X < 0.0 then
G := -G;
end if;
return Float_Type (G);
end Arcsin;
function Arcsin (X, Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Arcsin (X) * Cycle * One_Over_Two_Pi;
end Arcsin;
function Arccos (X : Float_Type) return Float_Type is
A : Float := Float (X);
Y : Float := Float (abs X);
F : Float;
G : Float;
I : Integer;
begin
if Y > 1.0 then
raise Argument_Error;
end if;
----Very small arguments are equal to the answer.
if Y < Low_Threshold then
I := 1;
G := Y;
else
----Compute the intermediate G value differently depending upon the
-- magnitude of the argument.
if Y > 0.5 then
I := 0;
F := ((0.5 - Y) + 0.5) / 2.0;
Y := Float (Sqrt (Float_Type (F)));
Y := -(Y + Y);
G := F;
else
I := 1;
G := Y * Y;
end if;
----Approximate our result.
F := ((((-0.69674_57344_73506_46411E+00 * G +
0.10152_52223_38064_63645E+02) * G -
0.39688_86299_75048_77339E+02) * G +
0.57208_22787_78917_31407E+02) * G -
0.27368_49452_41642_55994E+02) * G
/ (((((G - 0.23823_85915_36702_38830E+02) * G +
0.15095_27084_10306_04719E+03) * G -
0.38186_30336_17501_49284E+03) * G +
0.41714_43024_82604_12556E+03) * G -
0.16421_09671_44985_60795E+03);
G := Y + Y * F;
end if;
----Adjust our answer to a different quadrant if necessary.
if A < 0.0 then
if I = 1 then
G := 8#0.62207_73250_42055_06043# +
(8#0.62207_73250_42055_06043# + G);
else
G := 8#1.44417_66521_04132_14107# +
(8#1.44417_66521_04132_14107# + G);
end if;
else
if I = 1 then
G := 8#0.62207_73250_42055_06043# +
(8#0.62207_73250_42055_06043# - G);
else
G := -G;
end if;
end if;
----Return the result.
return Float_Type (G);
end Arccos;
function Arccos (X, Cycle : Float_Type) return Float_Type is
begin
if (abs X > 1.0) or (Cycle <= 0.0) then
raise Argument_Error;
end if;
return Arccos (X) * Cycle * One_Over_Two_Pi;
end Arccos;
function Arctan (Y : Float_Type;
X : Float_Type := 1.0) return Float_Type is
A : Float;
F : Float;
G : Float;
R : Float;
N : Natural range 0 .. 3;
begin
if X = 0.0 then
if Y = 0.0 then
raise Argument_Error;
elsif Y > 0.0 then
return Pi_Over_2;
else
return -Pi_Over_2;
end if;
end if;
if Y = 0.0 then
if X > 0.0 then
return 0.0;
else
return Pi;
end if;
end if;
----Reduce the argument to a small value close to zero.
if X = 1.0 then
A := Float (Y);
else
A := Float (Y) / Float (X);
end if;
F := abs A;
if F > 1.0 then
F := 1.0 / F;
N := 2;
else
N := 0;
end if;
if F > 8#0.21114_12136_47546_52614# then
F := (((8#0.56663_65641_30231_25163_55# * F - 0.5) - 0.5) + F) /
(8#1.56663_65641_30231_25163_55# + F);
N := N + 1;
end if;
----Approximate the result.
if abs F >= Low_Threshold then
G := F * F;
R := (((-0.83758_29936_81500_59274E+00 * G -
0.84946_24035_13206_83534E+01) * G -
0.20505_85519_58616_51981E+02) * G -
0.13688_76889_41919_26929E+02) * G
/ ((((G + 0.15024_00116_00285_76121E+02) * G +
0.59578_43614_25973_44465E+02) * G +
0.86157_34959_71302_42515E+02) * G +
0.41066_30668_25757_81263E+02);
F := F + F * R;
end if;
----Correct the result based upon the initial range reduction performed above.
if N > 1 then
F := -F;
end if;
case N is
when 0 =>
null;
when 1 =>
F := F + 0.52359_87755_98298_87308;
when 2 =>
F := F + 1.57079_63267_94896_61923;
when 3 =>
F := F + 1.04719_75511_96597_74615;
end case;
----Return a result of the correct sign.
if A < 0.0 then
F := -F;
end if;
return Float_Type (F);
end Arctan;
function Arctan (Y : Float_Type;
X : Float_Type := 1.0;
Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Arctan (Y, X) * Cycle * One_Over_Two_Pi;
end Arctan;
function Arccot (X : Float_Type;
Y : Float_Type := 1.0) return Float_Type is
begin
return Arctan (Y => Y, X => X);
end Arccot;
function Arccot (X : Float_Type;
Y : Float_Type := 1.0;
Cycle : Float_Type) return Float_Type is
begin
if Cycle <= 0.0 then
raise Argument_Error;
end if;
return Arccot (X, Y) * Cycle * One_Over_Two_Pi;
end Arccot;
function Sinh (X : Float_Type) return Float_Type is
A : Float := Float (X);
Y : Float := abs A;
W : Float;
Z : Float;
begin
----SINH for "large" arguments is computed this way.
if Y > 1.0 then
----Base-2 machines can compute things the easy way for "small" arguments.
if Y <= Log_Float_Large then
Z := Float (Exp (Float_Type (Y)));
Z := (Z * 0.5 - 0.5 / Z);
if A < 0.0 then
Z := -Z;
end if;
return Float_Type (Z);
end if;
----Reduce the argument.
W := Y - 8#0.542714#;
if W > Log_Float_Large then
raise Constraint_Error;
end if;
----Now take the EXP of the argument and compute the result with corrections
-- based on the machine radix.
Z := Float (Exp (Float_Type (W)));
Z := Z + 0.13830_27787_96019_02638E-04 * Z;
----Return the result with the correct sign.
if A < 0.0 then
Z := -Z;
end if;
return Float_Type (Z);
----Compute the result for "small" arguments here.
else
if Y < Low_Threshold then
return Float_Type (A);
end if;
Y := A * A;
W := (((-0.78966_12741_73570_99479E+00 * Y -
0.16375_79820_26307_51372E+03) * Y -
0.11563_52119_68517_68270E+05) * Y -
0.35181_28343_01771_17881E+06) * Y
/ (((Y - 0.27773_52311_96507_01667E+03) * Y +
0.36162_72310_94218_36460E+05) * Y -
0.21108_77005_81062_71242E+07);
W := A + A * W;
return Float_Type (W);
end if;
end Sinh;
function Cosh (X : Float_Type) return Float_Type is
A : Float := Float (X);
Y : Float := abs A;
W : Float;
Z : Float;
begin
----Base-2 machines can compute things the easy way for "small" arguments.
if Y <= Log_Float_Large then
Z := Float (Exp (Float_Type (Y)));
Z := (Z * 0.5 + 0.5 / Z);
return Float_Type (Z);
end if;
----Reduce the argument.
W := Y - 8#0.542714#;
if W > Log_Float_Large then
raise Constraint_Error;
end if;
----Now take the EXP of the argument and compute the result
Z := Float (Exp (Float_Type (W)));
Z := Z + 0.13830_27787_96019_02638E-04 * Z;
----Return the result.
return Float_Type (Z);
end Cosh;
function Tanh (X : Float_Type) return Float_Type is
A : Float := Float (X);
F : Float := abs A;
G : Float;
Result : Float;
Tanh_Xbig : constant Float := 1.90615474653985E+01;
-- Tanh_Xbig :=
-- (Ln_2 +
-- Float (Float'Machine_Mantissa + 1) * Ln (Float (Float'Machine_Radix)))
-- / 2.0;
begin
----Bigger arguments tend toward 1.0.
if F > Tanh_Xbig then
Result := 1.0;
----Middling big arguments can be best approximated this way.
elsif F > 0.54930_61443_34054_84570 then
Result := 0.5 - 1.0 / Float (Exp (Float_Type (F + F)) + 1.0);
Result := Result + Result;
----Very small arguments result in themselves.
elsif F < Low_Threshold then
Result := F;
----All other arguments get approximated here.
else
G := F * F;
Result := ((-0.96437_49277_72254_69787E+00 * G -
0.99225_92967_22360_83313E+02) * G -
0.16134_11902_39962_28058E+04) * G
/ (((G + 0.11274_47438_05349_49335E+03) * G +
0.22337_72071_89623_12926E+04) * G +
0.48402_35707_19886_88686E+04);
Result := F + F * Result;
end if;
----Return a result with the correct sign.
if A < 0.0 then
Result := -Result;
end if;
return Float_Type (Result);
end Tanh;
function Coth (X : Float_Type) return Float_Type is
begin
return 1.0 / Tanh (X);
end Coth;
function Arcsinh (X : Float_Type) return Float_Type is
begin
return Log (X + Sqrt (X * X + 1.0));
end Arcsinh;
function Arccosh (X : Float_Type) return Float_Type is
begin
return Log (X + Sqrt (X * X - 1.0));
end Arccosh;
function Arctanh (X : Float_Type) return Float_Type is
begin
return Log ((1.0 + X) / (1.0 - X)) / 2.0;
end Arctanh;
function Arccoth (X : Float_Type) return Float_Type is
begin
return Arctanh (1.0 / X);
end Arccoth;
begin
----Set up arrays of values for the Expx routine. See pages 98 and 106 of
-- the William J. Cody book for details of what and why.
declare
type Natural_Array is array (Positive range <>) of Natural;
procedure Expx_Stuff (A1 : in out Float;
A2 : in out Float;
N1 : Natural_Array) is
----N1/2/3 represent a floating point value. There are 3*30=90
-- bits of fraction in there. We take those 90 bits and we
-- place the first Float'Machine_Mantissa worth of them into A1
-- and we stick "the rest" into A2. No way to do this at
-- compile time.
Bits : Integer;
Mask : Integer;
begin
Bits := Float'Machine_Mantissa;
A1 := Float (N1 (1)) / 2.0 ** 30;
Mask := 2 ** (60 - Bits);
A1 := A1 + Float (N1 (2) - N1 (2) rem Mask) / 2.0 ** 60;
A2 := Float (N1 (2) rem Mask) / 2.0 ** 60 + Float (N1 (3)) / 2.0 ** 90;
end Expx_Stuff;
begin
Expx_A1 (1) := 1.0;
Expx_A1 (3) := 8#0.72540_30671_75644_41622_73201_32727#;
Expx_A1 (5) := 8#0.65642_37462_55323_53257_20715_15057#;
Expx_A1 (7) := 8#0.61263_45204_25240_66655_64761_25503#;
Expx_A1 (9) := 8#0.55202_36314_77473_63110_21313_73047#;
Expx_A1 (11) := 8#0.51377_32652_33052_11616_50345_72776#;
Expx_A1 (13) := 8#0.46033_76024_30667_05226_75065_32214#;
Expx_A1 (15) := 8#0.42712_70170_76521_36556_33737_10612#;
Expx_A1 (17) := 8#0.4#;
Expx_Stuff (Expx_A1 (2), Expx_A2 (1),
(8#75222_57505#, 8#22220_66302#, 8#61734_72062#));
Expx_Stuff (Expx_A1 (4), Expx_A2 (2),
(8#70146_33673#, 8#02522_47021#, 8#04062_61124#));
Expx_Stuff (Expx_A1 (6), Expx_A2 (3),
(8#63422_21405#, 8#21760_44016#, 8#17421_53016#));
Expx_Stuff (Expx_A1 (8), Expx_A2 (4),
(8#57204_24347#, 8#65401_41553#, 8#72504_02177#));
Expx_Stuff (Expx_A1 (10), Expx_A2 (5),
(8#53254_07672#, 8#44124_12254#, 8#31114_01243#));
Expx_Stuff (Expx_A1 (12), Expx_A2 (6),
(8#47572_46230#, 8#11064_10432#, 8#66404_42174#));
Expx_Stuff (Expx_A1 (14), Expx_A2 (7),
(8#44341_72334#, 8#72542_16063#, 8#30176_55544#));
Expx_Stuff (Expx_A1 (16), Expx_A2 (8),
(8#41325_30331#, 8#74611_03661#, 8#23056_22556#));
end;
end Generic_Elementary_Functions;
nblk1=37
nid=0
hdr6=6e
[0x00] rec0=14 rec1=00 rec2=01 rec3=02a
[0x01] rec0=1a rec1=00 rec2=02 rec3=04a
[0x02] rec0=03 rec1=00 rec2=37 rec3=01c
[0x03] rec0=22 rec1=00 rec2=03 rec3=03a
[0x04] rec0=01 rec1=00 rec2=36 rec3=012
[0x05] rec0=20 rec1=00 rec2=04 rec3=022
[0x06] rec0=00 rec1=00 rec2=35 rec3=006
[0x07] rec0=1b rec1=00 rec2=05 rec3=034
[0x08] rec0=2b rec1=00 rec2=06 rec3=05e
[0x09] rec0=01 rec1=00 rec2=34 rec3=000
[0x0a] rec0=22 rec1=00 rec2=07 rec3=02e
[0x0b] rec0=01 rec1=00 rec2=33 rec3=000
[0x0c] rec0=20 rec1=00 rec2=08 rec3=07e
[0x0d] rec0=04 rec1=00 rec2=32 rec3=00a
[0x0e] rec0=29 rec1=00 rec2=09 rec3=00a
[0x0f] rec0=01 rec1=00 rec2=31 rec3=002
[0x10] rec0=2c rec1=00 rec2=0a rec3=09c
[0x11] rec0=01 rec1=00 rec2=2f rec3=008
[0x12] rec0=21 rec1=00 rec2=30 rec3=02a
[0x13] rec0=02 rec1=00 rec2=0b rec3=00e
[0x14] rec0=22 rec1=00 rec2=0c rec3=014
[0x15] rec0=01 rec1=00 rec2=2e rec3=024
[0x16] rec0=26 rec1=00 rec2=0d rec3=016
[0x17] rec0=00 rec1=00 rec2=2d rec3=008
[0x18] rec0=2a rec1=00 rec2=0e rec3=028
[0x19] rec0=02 rec1=00 rec2=2c rec3=000
[0x1a] rec0=21 rec1=00 rec2=0f rec3=01a
[0x1b] rec0=00 rec1=00 rec2=2b rec3=008
[0x1c] rec0=2c rec1=00 rec2=10 rec3=034
[0x1d] rec0=02 rec1=00 rec2=2a rec3=004
[0x1e] rec0=1f rec1=00 rec2=11 rec3=024
[0x1f] rec0=00 rec1=00 rec2=29 rec3=008
[0x20] rec0=29 rec1=00 rec2=12 rec3=010
[0x21] rec0=00 rec1=00 rec2=28 rec3=010
[0x22] rec0=21 rec1=00 rec2=13 rec3=020
[0x23] rec0=00 rec1=00 rec2=27 rec3=00c
[0x24] rec0=2d rec1=00 rec2=14 rec3=008
[0x25] rec0=1b rec1=00 rec2=15 rec3=072
[0x26] rec0=2e rec1=00 rec2=16 rec3=02a
[0x27] rec0=1b rec1=00 rec2=17 rec3=00e
[0x28] rec0=26 rec1=00 rec2=18 rec3=054
[0x29] rec0=29 rec1=00 rec2=19 rec3=04e
[0x2a] rec0=1f rec1=00 rec2=1a rec3=012
[0x2b] rec0=28 rec1=00 rec2=1b rec3=01e
[0x2c] rec0=25 rec1=00 rec2=1c rec3=03a
[0x2d] rec0=22 rec1=00 rec2=1d rec3=010
[0x2e] rec0=27 rec1=00 rec2=1e rec3=040
[0x2f] rec0=00 rec1=00 rec2=26 rec3=05a
[0x30] rec0=22 rec1=00 rec2=1f rec3=00a
[0x31] rec0=00 rec1=00 rec2=25 rec3=00a
[0x32] rec0=2d rec1=00 rec2=20 rec3=01e
[0x33] rec0=15 rec1=00 rec2=21 rec3=010
[0x34] rec0=00 rec1=00 rec2=24 rec3=004
[0x35] rec0=11 rec1=00 rec2=22 rec3=088
[0x36] rec0=09 rec1=00 rec2=23 rec3=000
tail 0x215016e8481c04de3c61d 0x42a00066462061e03