|
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: 110870 (0x1b116)
Types: TextFile
Names: »B«
└─⟦85b835f43⟧ Bits:30000549 8mm tape, Rational 1000, Xlib rev 6.00
└─ ⟦0c20f784e⟧ »DATA«
└─⟦1abbe589f⟧
└─⟦306851c02⟧
└─⟦this⟧
--/ if R1000 then
with Debug_Tools;
--/ elsif TeleGen2 and then Unix then
--// with Float_Text_Io;
--/ end if;
with Text_Io;
with Xlbt_Arithmetic;
use Xlbt_Arithmetic;
package body Trig is
------------------------------------------------------------------------------
-- Date - /Name/ Comment
--
-- 01-FEB-89 - /GEB/ Created based upon the book
-- - Software Manual for the Elementary Functions by
-- - William J. Cody, Jr. and William Waite
-- - Prentice-Hall Series in Computational Mathematics, 1980
-- 05-JUL-90 - /GEB/ Add structured comments, port to TeleGen2/TeleSoft.
------------------------------------------------------------------------------
--/ if not (TeleGen2 and then Unix) then
package Float_Text_Io is new Text_Io.Float_Io (Float);
--/ end if;
type Float_Array is array (Natural32 range <>) of Float;
Degrees_To_Radians : constant Float := Pi / 180.0;
Sqrt_One_Tenth : constant Float :=
0.316227766016837933199889354443271853371955513727833;
Sqrt_One_Half : constant Float :=
0.7071067811865475244008443621048490392848359375608;
Ln_Float_Large : Float;
Ln_Float_Small : Float;
Ln_2 : constant Float :=
0.693147180559945309417232121458176568075500133107812;
Log10_E : constant Float :=
0.434294481903251827651128918916605082294397005627727;
Ln_10 : constant Float :=
2.30258509299404568401799145468436420760110148794495;
Expx_A1 : Float_Array (1 .. 17);
Expx_A2 : Float_Array (1 .. 8);
Tanh_Xbig : Float;
Cot_Eps1 : Float;
--\f
function Truncate (X : Float) return Float is
------------------------------------------------------------------------------
-- Chop off any non-integral portion of X. We do this by converting it to
-- an "integer" form and then back again without rounding.
------------------------------------------------------------------------------
Result : Float;
Fr : Float := abs X;
Ex : Integer32;
Tmp : Integer32;
Delt : Integer32;
begin
----Anything less than 1.0 in magnitude is 0.0 for our result.
if Fr < 1.0 then
return 0.0;
end if;
----Anything < Integer32'Last can be done using Integer32s.
if Fr < Float (Integer32'Last) then
Tmp := Integer32 (Fr);
if Float (Tmp) > Fr then
Tmp := Tmp - 1;
end if;
if X >= 0.0 then
return Float (Tmp);
else
return -Float (Tmp);
end if;
end if;
----Separate the fraction and the exponent in the number. Each time we take
-- some bit out of the fraction (in order to put them into our result) we
-- will take out out Delt Float'Machine_Radix "digits" worth. This guarantees
-- that we don't overflow our Integer32 temporary.
Fraction_Exponent (abs X, Fr, Ex);
case Float'Machine_Radix is
when 2 =>
Delt := (Integer32'Size - 1);
when 4 =>
Delt := (Integer32'Size - 1) / 2;
when 8 =>
Delt := (Integer32'Size - 1) / 3;
when 16 =>
Delt := (Integer32'Size - 1) / 4;
when 10 =>
Delt := (Integer32'Size - 1) / 4;
when others =>
Delt := 1;
end case;
----Loop extracting "large" chunks of bits from the fraction. Keep placing
-- them into the result. We will be extracting only integral bits.
Result := 0.0;
while Ex > Delt loop
Ex := Ex - Delt;
Fr := Fr * Float (Float'Machine_Radix) ** Integer (Delt);
Tmp := Integer32 (Fr);
if Float (Tmp) > Fr then -- Reverse any rounding affects.
Tmp := Tmp - 1;
end if;
Fr := Fr - Float (Tmp);
Result := Result * Float (Float'Machine_Radix) ** Integer (Delt) +
Float (Tmp);
end loop;
----Remove the last few Ex Float'Machine_Radix "digits" from the fraction.
-- This removes the last of the integral portion of the original number.
Fr := Fr * Float (Float'Machine_Radix) ** Integer (Ex);
Tmp := Integer32 (Fr);
if Float (Tmp) > Fr then
Tmp := Tmp - 1; -- Reverse any rounding affects.
end if;
Result := Result * Float (Float'Machine_Radix) ** Integer (Ex) +
Float (Tmp);
----Return a result with the correct sign.
if X < 0.0 then
return -Result;
else
return Result;
end if;
end Truncate;
--\f
function Truncate (X : Float) return Integer32 is
------------------------------------------------------------------------------
-- Convert the Float to an Integer32; truncate the fraction; do not round.
------------------------------------------------------------------------------
I : Integer32 := Integer32 (X);
begin
if I >= 0 then
if Float (I) > X then
return I - 1;
else
return I;
end if;
else
if Float (I) < X then
return I + 1;
else
return I;
end if;
end if;
end Truncate;
--\f
function Round (X : Float) return Float is
------------------------------------------------------------------------------
-- Chop off any non-integral portion of X. We do this by converting it to
-- an "integer" form and then back again with rounding.
------------------------------------------------------------------------------
Result : Float;
Fr : Float := abs X;
Ex : Integer32;
Tmp : Integer32;
Delt : Integer32;
begin
----Anything < Integer32'Last can be done using Integer32s.
if Fr < Float (Integer32'Last) then
Tmp := Integer32 (Fr);
if X >= 0.0 then
return Float (Tmp);
else
return -Float (Tmp);
end if;
end if;
----Separate the fraction and the exponent in the number. Each time we take
-- some bit out of the fraction (in order to put them into our result) we
-- will take out out Delt Float'Machine_Radix "digits" worth. This guarantees
-- that we don't overflow our Integer32 temporary.
Fraction_Exponent (abs X, Fr, Ex);
case Float'Machine_Radix is
when 2 =>
Delt := (Integer32'Size - 1);
when 4 =>
Delt := (Integer32'Size - 1) / 2;
when 8 =>
Delt := (Integer32'Size - 1) / 3;
when 16 =>
Delt := (Integer32'Size - 1) / 4;
when 10 =>
Delt := (Integer32'Size - 1) / 4;
when others =>
Delt := 1;
end case;
----Loop extracting "large" chunks of bits from the fraction. Keep placing
-- them into the result. We will be extracting only integral bits.
Result := 0.0;
while Ex > Delt loop
Ex := Ex - Delt;
Fr := Fr * Float (Float'Machine_Radix) ** Integer (Delt);
Tmp := Integer32 (Fr);
if Float (Tmp) > Fr then -- Reverse any rounding affects.
Tmp := Tmp - 1;
end if;
Fr := Fr - Float (Tmp);
Result := Result * Float (Float'Machine_Radix) ** Integer (Delt) +
Float (Tmp);
end loop;
----Remove the last few Ex Float'Machine_Radix "digits" from the fraction.
-- This removes the last of the integral portion of the original number.
Fr := Fr * Float (Float'Machine_Radix) ** Integer (Ex);
Tmp := Integer32 (Fr); -- Allow rounding to occur.
Result := Result * Float (Float'Machine_Radix) ** Integer (Ex) +
Float (Tmp);
----Return a result with the correct sign.
if X < 0.0 then
return -Result;
else
return Result;
end if;
end Round;
--\f
function Round (X : Float) return Integer32 is
begin
return Integer32 (X);
end Round;
--\f
procedure Fraction_Exponent (X : Float;
Fraction : out Float;
Exponent : out Integer32) is
------------------------------------------------------------------------------
-- Separate the number into a fraction (0.5 <= abs F < 1.0) and an exponent
-- such that X == Fraction * Float'Machine_Radix ** Exponent
------------------------------------------------------------------------------
F : Float := abs X;
E : Integer32 := 0;
begin
----A zero is 0.0 with 0 exponent.
if X = 0.0 then
Fraction := 0.0;
Exponent := 0;
return;
end if;
----Reduce values > 1.0 to values < 1.0. We scale by Machine_Radix so there
-- is no roundoff error whatsoever.
if F >= 1.0 then
while F >= Float (Float'Machine_Radix) ** 50 loop
E := E + 50;
F := F / Float (Float'Machine_Radix) ** 50;
end loop;
while F >= Float (Float'Machine_Radix) ** 10 loop
E := E + 10;
F := F / Float (Float'Machine_Radix) ** 10;
end loop;
while F >= 1.0 loop
E := E + 1;
F := F / Float (Float'Machine_Radix);
end loop;
----Increase values < 1.0/Machine_Radix to values >= 1.0/Machine_Radix. We
-- scale by Machine_Radix so there is no roundoff error whatsoever.
elsif F < 1.0 / Float (Float'Machine_Radix) then
while F < 1.0 / Float (Float'Machine_Radix) ** 50 loop
E := E - 50;
F := F * Float (Float'Machine_Radix) ** 50;
end loop;
while F < 1.0 / Float (Float'Machine_Radix) ** 10 loop
E := E - 10;
F := F * Float (Float'Machine_Radix) ** 10;
end loop;
while F < 1.0 / Float (Float'Machine_Radix) loop
E := E - 1;
F := F * Float (Float'Machine_Radix);
end loop;
end if;
----Return our results; preserve signs.
if X < 0.0 then
Fraction := -F;
else
Fraction := F;
end if;
Exponent := E;
end Fraction_Exponent;
--\f
function Odd (X : Float) return Boolean is
------------------------------------------------------------------------------
-- Returns TRUE if Aint(X) is an ODD value. Any value < 1.0 is even
-- because it is less than 1;
-- and any value > Float'Machine_Radix ** Float'Machine_Mantissa is
-- even because the even/odd bit has been lost due to insufficient
-- precision in the Float representation.
------------------------------------------------------------------------------
Y : Float := abs X;
Z : Float;
Shift : Integer32 := 0;
Int : Integer32;
begin
----If the value is too small to have an Integer32 value or if it has lost the
-- 1-bit that determines odd-ness then we are de-facto even.
if Y < 1.0 or else
Y >= Float (Float'Machine_Radix) ** Float'Machine_Mantissa then
return False;
end if;
----If necessary we will shift the number around until we can Integer32() it.
if Y >= Float (Integer32'Last) then
----Determine our shift amount. Guarantee that Y mod Z is < Integer32'Last.
case Float'Machine_Radix is
when 2 =>
Z := 2.0 ** (Integer32'Size - 2);
when 4 =>
Z := 4.0 ** ((Integer32'Size - 2) / 2);
when 8 =>
Z := 8.0 ** ((Integer32'Size - 2) / 3);
when 16 =>
Z := 16.0 ** ((Integer32'Size - 1) / 4);
when 10 =>
Z := 10.0 ** 9;
when others =>
raise Program_Error;
end case;
----Keep shifting down until we are small enough to fit in an Integer32 when
-- we are Integer32()'ed. Keep a count of the number of shifts we perform.
while Y >= Float (Integer32'Last) loop
Shift := Shift + 1;
Y := Y / Z;
end loop;
----Now recover those shifts. Remove successive groups of bits until we are
-- left with just the bottom Integer32'Size (more or less) bits from the
-- original Float.
while Shift > 0 loop
Int := Integer32 (Y);
if Float (Int) > Y then
Int := Int - 1;
end if;
Y := Y - Float (Int);
Y := Y * Z;
end loop;
end if;
----Check the last few bits and see if they constitute odd'ness.
Int := Integer32 (Y);
if Float (Int) > Y then
Int := Int - 1;
end if;
if Int rem 2 /= 0 then
return True;
else
return False;
end if;
end Odd;
--\f
function Sqrt (A : Float) return Float is
------------------------------------------------------------------------------
-- Square-Root - returns the positive root of positive numbers
-- - raises Constraint_Error for negative number
-- Based upon elements from:
-- Computer Approximations by John F. Hart et. al.; Krieger Publishing 1968
-- pg 89-96
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 17-34
------------------------------------------------------------------------------
Adj : Float;
F : Float;
N : Integer32;
Y : Float;
Bits : Natural32;
begin
----We don't do negative numbers and sqrt(0.0) = 0.0.
if A <= 0.0 then
if A < 0.0 then
raise Constraint_Error;
end if;
return 0.0;
end if;
----Break A into a fraction and an exponent.
Fraction_Exponent (A, F, N);
----Compute an approximation of the answer for sqrt(f).
case Float'Machine_Radix is
when 2 =>
Y := 0.41731 + 0.59016 * F;
when others =>
Y := 0.580661 + F * 0.5 - 0.086462 / (F + 0.175241);
end case;
----Now we put back the exponential information.
if N rem 2 /= 0 then
case Float'Machine_Radix is
when 2 =>
Y := Y * 8#0.55202_36314_77473_63110#;
when 4 =>
Y := Y * 4#0.2#;
when 8 =>
Y := Y * 8#0.26501_17146_37635_71444#;
when 10 =>
Y := Y * 0.31622_77660_16837_93320;
when 16 =>
Y := Y * 16#0.4#;
when others =>
raise Program_Error;
end case;
N := N + 1;
end if;
N := N / 2;
Y := Y * Float (Float'Machine_Radix) ** Integer (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) **
Integer (N - Float'Machine_Mantissa);
if not (Float'Machine_Rounds or else Y = Y + Adj / 2.0) or else
not (not Float'Machine_Rounds or else
Y = Y + Adj / 2.00001) or else
not (Y /= Y + Adj) or else not (Y = (Y + Adj) - Adj) then
raise Program_Error;
end if;
----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 =>
----We got an agument that was at/close-to Float'Large
-- Reduce Y back by Adj and return it.
Y := Y - Adj;
return Y;
end;
while F > A loop
Y := Y - Adj;
F := Y * Y;
end loop;
end if;
----Return the correct result.
return Y;
end Sqrt;
--\f
function Cbrt (A : Float) return Float is
------------------------------------------------------------------------------
-- Cube-Root - returns the real root
--
-- Based upon elements from:
-- Computer Approximations by John F. Hart et. al.; Krieger Publishing 1968
-- pg 89-96
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 17-34
------------------------------------------------------------------------------
Aa : Float := abs A;
Adj : Float;
F : Float;
N : Integer32;
Y : Float;
K : Integer32;
Bits : Natural32;
begin
----cbrt(0.0) = 0.0.
if A = 0.0 then
return 0.0;
end if;
----Break A into a fraction and an exponent.
Fraction_Exponent (Aa, F, N);
----Compute an approximation of the answer.
Y := (((-0.82075_488E-01 * F) + 0.42038_11298) * F + 0.78999_824) * F +
0.43261_46E-01;
Y := Y / (F + 0.17178_3523);
K := N rem 3;
if K > 0 then
----1/cbrt(3) = .69336.....
K := 3 - K;
Y := Y * 7.937005259840997373758528196362E-001 ** Integer (K);
N := N + K;
elsif K < 0 then
Y := Y * 7.937005259840997373758528196362E-001 ** Integer (abs K);
N := N - K;
end if;
N := N / 3;
Y := Y * Float (Float'Machine_Radix) ** Integer (N);
----Each time we do the equation we double our correct digits. We need
-- Float_Digits worth correct. Loop until this is true.
Bits := 8;
while Bits < Float'Size loop
Y := (Y + Y + Aa / (Y * Y)) / 3.0;
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 * Y;
if F /= Aa then
Adj := Float (Float'Machine_Radix) **
Integer (N - Float'Machine_Mantissa);
if not (Float'Machine_Rounds or else Y = Y + Adj / 2.0) or else
not (not Float'Machine_Rounds or else
Y = Y + Adj / 2.00001) or else
not (Y /= Y + Adj) or else
not (Y = (Y + Adj) - Adj) then
raise Program_Error;
end if;
----Adj is so small that it only diddles the bottom bit of the
-- mantissa; we so assert.
begin
while F < Aa loop
Y := Y + Adj;
F := Y * Y * Y;
end loop;
exception
when Numeric_Error =>
----We got an agument that was at/close-to Float'Large
-- Reduce Y back by Adj and return it.
Y := Y - Adj;
return Y;
end;
while F > Aa loop
Y := Y - Adj;
F := Y * Y * Y;
end loop;
end if;
----Return the correct result.
if A < 0.0 then
return -Y;
else
return Y;
end if;
end Cbrt;
--\f
function Sin (A : Float) return Float is
R : Float;
F : Float;
G : Float;
Y : Float;
Xn : Float;
Tmp : Float;
Sign : Boolean;
begin
----Get Y := abs A where Rad is also in radians.
if In_Radians then
Y := A;
else
Y := A * Degrees_To_Radians;
end if;
if Y < 0.0 then
Y := -Y;
Sign := True;
else
Sign := False;
end if;
----Check for ridiculously large values of Y.
if Y >= Pi * Float (Float'Machine_Radix) **
(Float'Machine_Mantissa / 2) then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (Y * 0.31830_98861_83790_67154);
if Odd (Xn) then
Sign := not Sign;
end if;
----Compute the f portion of the argument.
Tmp := Truncate (Y);
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) then
F := ((Tmp - Xn * 8#3.1104#) + (Y - Tmp)) +
Xn * 8.9089_10206_76153_73566_17E-06;
else
F := ((Tmp - Xn * 8#3.11#) + (Y - Tmp)) - Xn * 9.6765_35897_93E-04;
end if;
----See if there is no point to computing sin(x).
if abs F < Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
if Sign then
return -F;
else
return F;
end if;
end if;
----Compute an answer.
G := F * F;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 51) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 51) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 51) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 51) or else
(Float'Machine_Radix = 10 and then Float'Machine_Mantissa >= 16) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 11) then
R := ((((((-0.73706_62775_07114_174E-012 * G +
0.16047_84463_23816_900E-09) * G -
0.25051_87088_34705_760E-07) * G +
0.27557_31642_12926_457E-05) * G -
0.19841_26982_32225_068E-03) * G +
0.8333_33333_27592_139E-02) * G -
0.16666_66666_66659_653E+00) * G;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 9) then
R := ((((-0.23868_34640_601E-07 * G +
0.27523_97106_775E-05) * G -
0.19840_83282_313E-03) * G +
0.83333_30720_556E-02) * G -
0.16666_66660_883E+00) * G;
else
R := (((0.26019_03036E-05 * G -
0.19807_41872E-03) * G +
0.83330_25139E-02) * G -
0.16666_65668E+00) * G;
end if;
----Now finish off the result.
R := F + F * R;
if Sign then
return -R;
else
return R;
end if;
end Sin;
--\f
function Cos (A : Float) return Float is
R : Float;
F : Float;
G : Float;
Y : Float;
Xn : Float;
Tmp : Float;
Sign : Boolean;
begin
----Get Y := abs A where Rad is also in radians.
if In_Radians then
Y := abs A;
else
Y := abs A * Degrees_To_Radians;
end if;
Y := Y + 1.57079_63267_94896_61923;
Sign := False;
----Check for ridiculously large values of Y.
if Y >= Pi * Float (Float'Machine_Radix) **
(Float'Machine_Mantissa / 2) then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (Y * 0.31830_98861_83790_67154);
if Odd (Xn) then
Sign := not Sign;
end if;
Xn := Xn - 0.5;
----Compute the f portion of the argument.
Tmp := Truncate (abs A);
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) then
F := ((Tmp - Xn * 8#3.1104#) + (abs A - Tmp)) +
Xn * 8.9089_10206_76153_73566_17E-06;
else
F := ((Tmp - Xn * 8#3.11#) + (abs A - Tmp)) -
Xn * 9.6765_35897_93E-04;
end if;
----See if there is no point to computing cos(x).
if abs F < Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
if Sign then
return -F;
else
return F;
end if;
end if;
----Compute an answer.
G := F * F;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 51) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 51) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 51) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 51) or else
(Float'Machine_Radix = 10 and then Float'Machine_Mantissa >= 16) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 11) then
R := ((((((-0.73706_62775_07114_174E-012 * G +
0.16047_84463_23816_900E-09) * G -
0.25051_87088_34705_760E-07) * G +
0.27557_31642_12926_457E-05) * G -
0.19841_26982_32225_068E-03) * G +
0.8333_33333_27592_139E-02) * G -
0.16666_66666_66659_653E+00) * G;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 9) then
R := ((((-0.23868_34640_601E-07 * G +
0.27523_97106_775E-05) * G -
0.19840_83282_313E-03) * G +
0.83333_30720_556E-02) * G -
0.16666_66660_883E+00) * G;
else
R := (((0.26019_03036E-05 * G -
0.19807_41872E-03) * G +
0.83330_25139E-02) * G -
0.16666_65668E+00) * G;
end if;
----Now finish off the result.
R := F + F * R;
if Sign then
return -R;
else
return R;
end if;
end Cos;
--\f
function Tan (A : Float) return Float is
Aa : Float;
F : Float;
X1 : Float;
Xnum : Float;
Xden : Float;
Even : Boolean;
Xn : Float;
begin
----Make sure we have the correct units.
if In_Radians then
Aa := A;
else
Aa := A * Degrees_To_Radians;
end if;
F := abs Aa;
----See if our argument is much too large to make any sense.
if F > Truncate (Float (Float'Machine_Radix) **
(Float'Machine_Mantissa / 2) * Pi_Over_2) then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (Aa * 0.63661_97723_67581_34308);
Even := not Odd (Xn);
----Compute f.
X1 := Truncate (Aa);
if Float'Machine_Radix /= 10 then
F := ((X1 - Xn * 8#1.4442#) + (Aa - X1)) +
Xn * 4.4544_55103_38706_86783_08E-06;
else
F := ((X1 - Xn * 1.57079) + (Aa - X1)) +
Xn * 6.32679489661923132169163975144E-06;
end if;
if not (abs F <= Pi_Over_4) then
raise Program_Error;
end if;
----See if f is very very small.
if abs F < Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
Xnum := F;
Xden := 1.0;
else
X1 := F * F;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 53) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 53) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 53) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 53) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 17) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 11) then
Xnum := (((-0.74836_34966_61206_5149E-05 * X1 +
0.28059_18241_16998_8906E-02) * X1 -
0.12828_34704_09574_3847E+00) * X1) * F + F;
Xden := ((-0.20844_80442_20387_0948E-03 * X1 +
0.23344_85282_20687_2802E-01) * X1 -
0.46161_68037_42904_8840E+00) * X1 +
0.5;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 8) then
Xnum := ((0.10751_54738_488E-02 * X1 -
0.11136_14403_566E+00) * X1) * F + F;
Xden := (0.15973_39213_300E-01 * X1 -
0.44469_47720_281E+00) * X1 +
0.5;
else
Xnum := (-0.95801_7723E-01 * X1) * F + F;
Xden := (0.97168_5835E-02 * X1 -
0.42913_5777E+00) * X1 +
0.5;
end if;
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 Xnum;
end Tan;
--\f
function Sec (A : Float) return Float is
begin
return 1.0 / Cos (A);
end Sec;
--\f
function Csc (A : Float) return Float is
begin
return 1.0 / Sin (A);
end Csc;
--\f
function Cot (A : Float) return Float is
Aa : Float;
F : Float;
X1 : Float;
Xnum : Float;
Xden : Float;
Even : Boolean;
Xn : Float;
begin
----Make sure we have the correct units.
if In_Radians then
Aa := A;
else
Aa := A * Degrees_To_Radians;
end if;
F := abs Aa;
----See if our argument is much too small to handle.
if F < Cot_Eps1 then
if Aa = 0.0 then
return Float'Last;
elsif Aa < 0.0 then
return Float'First;
else
return Float'Last;
end if;
end if;
----See if our argument is much too large to make any sense.
if F > Truncate (Float (Float'Machine_Radix) **
(Float'Machine_Mantissa / 2) * Pi_Over_2) then
return 0.0;
end if;
----Form XN and check for even/odd-ness.
Xn := Round (Aa * 0.63661_97723_67581_34308);
Even := not Odd (Xn);
----Compute f.
X1 := Truncate (Aa);
if Float'Machine_Radix /= 10 then
F := ((X1 - Xn * 8#1.4442#) + (Aa - X1)) +
Xn * 4.4544_55103_38706_86783_08E-06;
else
F := ((X1 - Xn * 1.57079) + (Aa - X1)) +
Xn * 6.32679489661923132169163975144E-06;
end if;
----See if f is very very small.
if abs F < Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
Xnum := F;
Xden := 1.0;
else
X1 := F * F;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 53) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 53) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 53) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 53) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 17) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 11) then
Xnum := (((-0.74836_34966_61206_5149E-05 * X1 +
0.28059_18241_16998_8906E-02) * X1 -
0.12828_34704_09574_3847E+00) * X1) * F + F;
Xden := ((-0.20844_80442_20387_0948E-03 * X1 +
0.23344_85282_20687_2802E-01) * X1 -
0.46161_68037_42904_8840E+00) * X1 +
0.5;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 8) then
Xnum := ((0.10751_54738_488E-02 * X1 -
0.11136_14403_566E+00) * X1) * F + F;
Xden := (0.15973_39213_300E-01 * X1 -
0.44469_47720_281E+00) * X1 +
0.5;
else
Xnum := (-0.95801_7723E-01 * X1) * F + F;
Xden := (0.97168_5835E-02 * X1 -
0.42913_5777E+00) * X1 +
0.5;
end if;
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 Xnum;
end Cot;
--\f
function Asin (A : Float) return Float is
Y : Float := abs A;
F : Float;
G : Float;
I : Integer32;
begin
----Very small arguments are equal to the answer.
if Y < Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa / 2) 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;
if Y > 1.0 then
raise Constraint_Error;
end if;
F := ((0.5 - Y) + 0.5) / 2.0;
Y := Sqrt (F);
Y := -(Y + Y);
G := F;
else
I := 0;
G := Y * Y;
end if;
----Approximate our result.
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 49) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 49) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 49) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 49) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 15) then
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);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 37) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 37) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 37) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 37) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 12) then
F := (((-0.65404_06899_93350_09E+00 * G +
0.59683_15761_77515_34E+01) * G -
0.13428_70791_34253_12E+02) * G +
0.85372_16436_67719_50E+01) * G
/ ((((G - 0.16429_55755_74951_70E+02) * G +
0.68729_59765_38088_06E+02) * G -
0.10362_27318_64014_80E+03) * G +
0.51223_29862_01096_91E+02);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 8) then
F := ((-0.59450_14419_3246E+00 * G +
0.29058_76237_4859E+01) * G -
0.27516_55529_0596E+01) * G
/ (((G - 0.10333_86707_2113E+02) * G +
0.24864_72896_9164E+02) * G -
0.16509_93320_2424E+02);
else
F := (-0.50440_0557E+00 * G +
0.93393_5835E+00) * G
/ ((G - 0.55484_6723E+01) * G +
0.56036_3004E+01);
end if;
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 not In_Radians then
G := G / Degrees_To_Radians;
end if;
if A < 0.0 then
return -G;
else
return G;
end if;
end Asin;
--\f
function Acos (A : Float) return Float is
Y : Float := abs A;
F : Float;
G : Float;
I : Integer32;
begin
----Very small arguments are equal to the answer.
if Y < Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa / 2) 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;
if Y > 1.0 then
raise Constraint_Error;
end if;
F := ((0.5 - Y) + 0.5) / 2.0;
Y := Sqrt (F);
Y := -(Y + Y);
G := F;
else
I := 1;
G := Y * Y;
end if;
----Approximate our result.
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 49) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 49) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 49) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 49) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 15) then
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);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 37) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 37) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 37) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 37) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 12) then
F := (((-0.65404_06899_93350_09E+00 * G +
0.59683_15761_77515_34E+01) * G -
0.13428_70791_34253_12E+02) * G +
0.85372_16436_67719_50E+01) * G
/ ((((G - 0.16429_55755_74951_70E+02) * G +
0.68729_59765_38088_06E+02) * G -
0.10362_27318_64014_80E+03) * G +
0.51223_29862_01096_91E+02);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 8) then
F := ((-0.59450_14419_3246E+00 * G +
0.29058_76237_4859E+01) * G -
0.27516_55529_0596E+01) * G
/ (((G - 0.10333_86707_2113E+02) * G +
0.24864_72896_9164E+02) * G -
0.16509_93320_2424E+02);
else
F := (-0.50440_0557E+00 * G +
0.93393_5835E+00) * G
/ ((G - 0.55484_6723E+01) * G +
0.56036_3004E+01);
end if;
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.
if not In_Radians then
G := G / Degrees_To_Radians;
end if;
return G;
end Acos;
--\f
function Atan (A : Float) return Float is
F : Float;
G : Float;
R : Float;
N : Natural32 range 0 .. 3;
begin
----Reduce the argument to a small value close to zero.
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 >= Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
G := F * F;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 51) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 51) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 51) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 51) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 16) then
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);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 11) then
R := ((-0.79439_12954_08336_251E+00 * G -
0.42744_49853_67930_329E+01) * G -
0.42743_26720_26241_096E+01) * G
/ (((G + 0.91978_93648_35039_806E+01) * G +
0.20517_13765_64218_456E+02) * G +
0.12822_98016_07919_841E+02);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 9) then
R := (-0.72002_68488_98E+00 * G -
0.14400_83448_74E+01) * G
/ ((G + 0.47522_25845_99E+01) * G +
0.43202_50389_19E+01);
else
R := (-0.50909_58253E-01 * G -
0.47083_25141E+00) * G
/ (G + 0.14125_00740E+01);
end if;
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 not In_Radians then
F := F / Degrees_To_Radians;
end if;
if A < 0.0 then
F := -F;
end if;
return F;
end Atan;
--\f
function Asec (A : Float) return Float is
begin
return Acos (1.0 / A);
end Asec;
--\f
function Acsc (A : Float) return Float is
begin
return Asin (1.0 / A);
end Acsc;
--\f
function Acot (A : Float) return Float is
begin
return Atan (1.0 / A);
end Acot;
--\f
function Ln (A : Float) return Float is
------------------------------------------------------------------------------
-- Ln - returns the natural logarithm of its positive argument
--
-- Based upon elements from
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 35-48
------------------------------------------------------------------------------
F : Float;
N : Integer32;
Znum : Float;
Zden : Float;
begin
----Base=10 systems use the log10 function.
if Float'Machine_Radix = 10 then
return Log10 (A) / Log10_E;
else
----Negative arguments are errors. A 0.0 argument has an "easy" answer.
if A < 0.0 then
raise Constraint_Error;
end if;
if A = 0.0 then
raise Numeric_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.
Fraction_Exponent (A, F, N);
case Float'Machine_Radix is
when 2 =>
null;
when 4 =>
N := N * 2;
while F < 0.5 loop
F := F + F;
N := N - 1;
end loop;
when 8 =>
N := N * 3;
while F < 0.5 loop
F := F + F;
N := N - 1;
end loop;
when 16 =>
N := N * 4;
while F < 0.5 loop
F := F + F;
N := N - 1;
end loop;
when others =>
raise Program_Error;
end case;
----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;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 49) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 49) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 49) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 49) then
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);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 33) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 33) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 33) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 33) then
Zden :=
((F * 0.44445_51510_98033_23E-02 -
0.63260_86623_38596_65E+00)
* F + 0.37339_16896_31608_66E+01) /
((F - 0.14312_35435_58853_24E+02)
* F + 0.44807_00275_57364_36E+02);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa >= 25) then
Zden := (F * 0.13600_95468_621E-01 - 0.46490_62303_464E+00) /
(F - 0.55788_73750_242E+01);
else
Zden := (-0.555270_74855E+00) / (F - 0.66327_18214E+01);
end if;
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 F;
end if;
end Ln;
--\f
function Log10 (A : Float) return Float is
------------------------------------------------------------------------------
-- Log10 - returns the base-10 logarithm of its positive argument
--
-- Based upon elements from
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 35-48
------------------------------------------------------------------------------
F : Float;
N : Integer32;
S : Float;
S2 : Float;
begin
----Base/=10 systems use the ln function.
if Float'Machine_Radix /= 10 then
return Ln (A) / Ln_10;
else
----Negative arguments are errors. A 0.0 argument has an "easy" answer.
if A < 0.0 then
raise Constraint_Error;
end if;
if A = 0.0 then
raise Numeric_Error;
end if;
----Extract the fraction and the exponent from the argument.
Fraction_Exponent (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_Tenth then
null;
else
N := N - 1;
F := F * 10.0;
end if;
----Now we evaluate a rational approximation of exp(z).
S := ((F - 0.5) - 0.5) / (F + 1.0);
S2 := S * S;
if Float'Machine_Mantissa >= 16 then
F := ((((-0.74101_07841_61919_23924E+00 * S2 +
0.10333_85715_14793_86456E+02)
* S2 - 0.39273_74102_03156_25018E+02)
* S2 + 0.55408_59120_41205_93097E+02)
* S2 - 0.26044_70024_05557_63612E+02)
/ (((((S2 - 0.19373_23458_32854_78594E+2)
* S2 + 0.10710_97891_15668_00857E+03)
* S2 - 0.24430_30353_41829_54205E+03)
* S2 + 0.24534_76188_68489_34790E+3)
* S2 - 0.89955_20778_81033_11704E+02);
elsif Float'Machine_Mantissa >= 11 then
F := (((-0.71433_38215_32264_273E+00 * S2 +
0.62503_65112_79083_731E+01)
* S2 - 0.13682_37024_15026_896E+02)
* S2 + 0.85167_31987_23885_403E+01)
/ ((((S2 - 0.13210_47835_01562_817E+02)
* S2 + 0.47925_25604_38733_968E+02)
* S2 - 0.64906_68274_09428_483E+02)
* S2 + 0.29415_75017_23226_173E+02);
elsif Float'Machine_Mantissa >= 8 then
F := ((-0.67358_16014_777E+00 * S2 + 0.31630_34915_570E+01)
* S2 - 0.29156_81437_901E+01)
/ (((S2 - 0.81908_00454_670E+01)
* S2 + 0.16966_98140_210E+02)
* S2 - 0.10070_40695_423E+02);
else
F := (-0.60368_24627E+00 * S2 + 0.10756_13712E+01) /
((S2 - 0.43144_78001E+01) * S2 + 0.37150_53570E+01);
end if;
S2 := S2 * F;
F := S * (Log10_E + S2);
----Now compute the result.
F := F + Float (N);
return F;
end if;
end Log10;
--\f
function Logx (X : Float; Y : Float) return Float is
begin
if Float'Machine_Radix /= 10 then
return Ln (Y) / Ln (X);
else
return Log10 (Y) / Log10 (X);
end if;
end Logx;
--\f
function Exp (A : Float) return Float is
------------------------------------------------------------------------------
-- Exp - returns E to the power represented by its argument
--
-- Based upon elements from Computer Approximations by John F. Hart et. al.
-- Krieger Publishing 1968
-- pg 97-105
------------------------------------------------------------------------------
N : Integer32;
Xn : Float;
G : Float;
Xi : Integer32;
X1 : Float;
X2 : Float;
Pz : Float;
Qz : Float;
begin
----Make sure that it is possible to represent the result. (too big/small)
if A > Ln_Float_Large then
raise Numeric_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;
----Handle non-base-10 machines here.
if Float'Machine_Radix /= 10 then
----Perform agrument range reduction.
N := Integer32
(A * 1.44269504088896340735992468100189213742664595602385);
Xn := Float (N);
Xi := Integer32 (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;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 57) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 57) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 57) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 57) then
Pz := (0.31555_19276_56846_46356E-04 * X1 +
0.75753_18015_94227_76666E-02) * X1 + 0.25;
Qz := ((0.75104_02839_98700_46114E-06 * X1 +
0.63121_89437_43985_03557E-03) * X1 +
0.56817_30269_85512_21787E-01) * X1 + 0.5;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 43) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 43) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 43) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 43) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 30) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 30) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 30) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 30) then
Pz := 0.59504_25497_7591E-02 * X1 + 0.24999_99999_9992E+00;
Qz :=
(0.29729_36368_2238E-03 * X1 + 0.53567_51764_522E-01) * X1 +
0.5;
else
Pz := 0.41602_88626_8E-02 * X1 + 0.24999_99995_0E+00;
Qz := 0.49987_17877_8E-01 * X1 + 0.5;
end if;
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 ** Integer (Xi);
G := G * 2.0 ** Integer (N);
return G;
----Handle base-10 machines here.
else
----Perform agrument range reduction.
N := Integer32
(A * 0.868588963806503655302257837833210164588794011255454);
Xn := Float (N);
Xi := Integer32 (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 * 1.151E+00) + X2) - Xn * 2.9254_64970_22842_009E-04;
----Now approximate exp(g)
X1 := G * G;
if Float'Machine_Mantissa >= 15 then
Pz := (0.42041_42681_37450_31524E+02 * X2 +
0.10097_41487_24273_91798E+05) * X2 +
0.33326_70292_26801_61137E+06;
Qz := ((X2 + 0.84124_35845_14154_54514E+03) * X2 +
0.75739_33461_59883_44395E+05) * X2 +
0.66653_40584_53603_22323E+06;
elsif Float'Machine_Mantissa >= 11 then
Pz := (0.33287_36465_16410_3E-01 * X1 +
0.14008_29975_62819_6E+02) * X1 +
0.50446_48895_05869_7E+03;
Qz := (X1 + 0.11209_40810_96616_9E+03) * X1 +
0.50445_48895_05869_7E+03;
elsif Float'Machine_Mantissa >= 8 then
Pz := 0.20041_43275_526E+02 * X1 + 0.84249_03867_900E+03;
Qz := (X2 + 0.18049_79288_462E+03) * X2 + 0.16849_80773_608E+04;
else
Pz := 0.83046_5413E-01 * X1 + 0.50034_9857E+01;
Qz := X1 + 0.10006_9975E+02;
end if;
Pz := Pz * G;
G := 0.5 + Pz / (Qz - Pz);
G := G + G;
----Compute the final result.
if N rem 2 /= 0 then
if N > 0 then
G :=
G * 3.16227766016837933199889354443271853371955514008093;
else
G :=
G / 3.16227766016837933199889354443271853371955514008093;
end if;
end if;
N := N / 2;
Xi := N / 2;
N := N - Xi;
G := G * 2.0 ** Integer (Xi);
G := G * 2.0 ** Integer (N);
return G;
end if;
end Exp;
--\f
function Exp10 (A : Float) return Float is
begin
return Expx (10.0, A);
end Exp10;
--\f
function Expx (X : Float; Y : Float) return Float is
M : Integer32;
F : Float;
R : Integer32;
P : Integer32;
Z : Float;
Rv : Float;
U1 : Float;
U2 : Float;
Y1 : Float;
Y2 : Float;
Iw1 : Integer32;
W : Float;
W1 : Float;
W2 : Float;
I : Integer32;
N : Integer32;
begin
----We use the conventions where (-f)**g is illegal and f**0.0 is 1.0.
if Y = 0.0 then
return 1.0;
elsif X < 0.0 then
raise Constraint_Error;
end if;
----Here we handle non-base-10 floating point machines.
if Float'Machine_Radix /= 10 then
Fraction_Exponent (X, F, M);
case Float'Machine_Radix is
when 2 =>
null;
when 4 =>
if F < 0.5 then
R := 1;
F := F + F;
else
R := 0;
end if;
when 8 =>
if F < 0.25 then
R := 2;
F := F + F;
F := F + F;
elsif F < 0.5 then
R := 1;
F := F + F;
else
R := 0;
end if;
when 16 =>
if F < 0.125 then
R := 3;
F := F + F;
F := F + F;
F := F + F;
elsif F < 0.25 then
R := 2;
F := F + F;
F := F + F;
elsif F < 0.5 then
R := 1;
F := F + F;
else
R := 0;
end if;
when others =>
raise Program_Error;
end case;
----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.
if Float'Machine_Radix /= 2 then
Z := F - Expx_A1 (P + 1);
Z := Z - Expx_A2 ((P + 1) / 2);
Z := Z / (F * 0.5 + Expx_A1 (P + 1) * 0.5);
else
Z := F - Expx_A1 (P + 1);
Z := Z - Expx_A2 ((P + 1) / 2);
Z := Z / (F + Expx_A1 (P + 1));
Z := Z + Z;
end if;
----Now compute the z*P(z**2) approximation.
F := Z * Z;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 51) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 51) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 51) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 51) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 37) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 37) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 37) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 37) then
Rv :=
(0.22338_24352_81541_8E-02 * F + 0.12499_99796_50060_8E-01) *
F +
0.83333_33333_41213_6E-01;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) then
Rv := 0.12506_48500_52E-01 * F + 0.83333_32862_45E-01;
else
Rv := 0.83357_541E-01;
end if;
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.
case Float'Machine_Radix is
when 2 =>
U1 := Float (M * 16 - P) * 0.0625;
when 4 =>
U1 := Float ((2 * M - R) * 16 - P) * 0.0625;
when 8 =>
U1 := Float ((3 * M - R) * 16 - P) * 0.0625;
when 16 =>
U1 := Float ((4 * M - R) * 16 - P) * 0.0625;
when others =>
raise Program_Error;
end case;
Y1 := Truncate (Y * 16.0) * 0.0625;
Y2 := Y - Y1;
W := U2 * 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 := 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.
case Float'Machine_Radix is
when 2 =>
M := Iw1 / 16 + I;
P := 16 * M - Iw1;
when 4 =>
N := Iw1 / 16 + I;
P := 16 * N - Iw1;
M := N / 2 + I;
R := 2 * M - N;
when 8 =>
N := Iw1 / 16 + I;
P := 16 * N - Iw1;
M := N / 3 + I;
R := 3 * M - N;
when 16 =>
N := Iw1 / 16 + I;
P := 16 * N - Iw1;
M := N / 4 + I;
R := 4 * M - N;
when others =>
raise Program_Error;
end case;
----Now do Z.
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 53) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 53) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 53) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 53) then
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;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 46) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 46) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 46) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 46) then
Z := (((((0.15077_40617_88142_382E-03 * W2 +
0.13330_81011_34082_075E-02) * W2 +
0.96181_17691_38724_104E-02) * W2 +
0.55504_10842_47568_661E-01) * W2 +
0.24022_65069_56777_522E+00) * W2 +
0.69314_71805_59937_815E+00) * W2;
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) then
Z := ((((0.13052_55159_42810E-02 * W2 +
0.96162_06595_83789E-02) * W2 +
0.55504_04881_30765E-01) * W2 +
0.24022_65061_44710E+00) * W2 +
0.69314_71805_56341E+00) * W2;
else
Z := ((0.54360_383E-01 * W2 +
0.24018_510E+00) * W2 +
0.69314_675E+00) * W2;
end if;
----Nearly there; one more step.
Z := Expx_A1 (P + 1) + Expx_A1 (P + 1) * Z;
if Float'Machine_Radix /= 2 then
Z := Z * 2.0 ** Integer (R);
end if;
I := M / 2;
M := M - I;
Z := Z * Float (Float'Machine_Radix) ** Integer (I);
Z := Z * Float (Float'Machine_Radix) ** Integer (M);
return Z;
----Here we handle base-10 floating point machines.
else
Fraction_Exponent (X, F, M);
----Compute the P index.
P := 1;
while F <= Expx_A1 (P + 2) loop
P := P + 2;
end loop;
----Compute z.
Z := F - Expx_A1 (P + 1);
Z := Z - Expx_A2 ((P + 1) / 2);
Z := Z / (F + Expx_A1 (P + 1));
----Now compute the z*P(z**2) approximation.
F := Z * Z;
if Float'Machine_Mantissa >= 15 then
Rv := ((-0.68650_74087_13682_76474E+00 * F +
0.35829_52156_67236_56966E+01) * F -
0.35184_51552_93711_75192E+01) * F
/ (((F - 0.89628_58940_90517_72560E+01) * F +
0.19666_45902_38905_35938E+02) * F -
0.12152_30114_43221_36016E+02);
elsif Float'Machine_Mantissa >= 9 then
Rv := (-0.62077_57571_87918_2E+00 * F +
0.12095_10556_17031_7E+01) * F
/ ((F - 0.46505_84419_61834_9E+01) * F +
0.41775_01464_74140_2E+01);
else
Rv := -0.48144_78516E+00 * F / (F - 0.16628_73334E+01);
end if;
U2 := Z * (Rv + 0.86858_89638_06503_65530);
----Now create W, our partial answer.
U1 := Float (M * 16 - P) * 0.0625;
Y1 := Truncate (Y * 10.0) * 0.1;
Y2 := Y - Y1;
W := U2 * Y + U1 * Y2;
W1 := Truncate (W * 10.0) * 0.1;
W2 := W - W1;
W := W1 + U1 * Y1;
W1 := Truncate (W * 10.0) * 0.1;
W2 := W2 + (W - W1);
W := Truncate (W2 * 10.0) * 0.1;
Iw1 := Truncate (10.0 * (W1 + W));
W2 := W2 - W;
if Iw1 > 10 * Float'Machine_Emax - 1 then
raise Numeric_Error; -- Result far too large
end if;
if Iw1 < 10 * Float'Machine_Emin + 1 then
return 0.0; -- Result far too small
end if;
if W2 > 0.0 then
Iw1 := Iw1 + 1;
W2 := W2 - 0.1;
end if;
if Iw1 < 0 then
I := 0;
else
I := 1;
end if;
----Compute the new p,m values.
M := Iw1 / 10 + I;
P := 10 * M - Iw1;
----Now do Z.
if Float'Machine_Mantissa >= 15 then
Z := ((((0.41071_57590_70777_89734E+0 * W2 +
0.35616_93889_90563_55495E+01) * W2 +
0.34363_74379_23829_31895E+02) * W2 +
0.28195_17958_42955_14217E+02) * W2 +
0.23539_70807_49261_19272E+03) * W2
/ ((((W2 - 0.10224_89787_60959_10327E+02) * W2 +
0.45994_87860_33907_54883E+02) * W2 -
0.10545_35294_64899_79088E+03) * W2 +
0.10223_16532_25538_32434E+03);
elsif Float'Machine_Mantissa >= 11 then
Z := (((-0.51359_30940_15587_9E+0 * W2 -
0.46710_07431_07040_0E+01) * W2 -
0.57975_91698_28398_6E+01) * W2 -
0.37713_99188_31947_3E+02) * W2
/ (((W2 - 0.63664_36505_45209_5E+01) * W2 +
0.16339_13385_89750_1E+02) * W2 -
0.16378_97856_54144_4E+02);
elsif Float'Machine_Mantissa >= 8 then
Z := ((0.68536_85048_84E+00 * W2 +
0.17805_77775_84E+01) * W2 +
0.82998_14853_06E+01) * W2
/ ((W2 - 0.33766_11203_26E+01) * W2 +
0.36045_63799_45E+01);
else
Z := ((-0.78964_7014E+00 * W2 -
0.21742_4847E+01) * W2 -
0.38883_9832E+01) * W2
/ (W2 - 0.16887_0925E+01);
end if;
----Nearly there; one more step.
Z := Expx_A1 (P + 1) + Expx_A1 (P + 1) * Z;
I := M / 2;
M := M - I;
Z := Z * 10.0 ** Integer (I);
Z := Z * 10.0 ** Integer (M);
return Z;
end if;
end Expx;
--\f
function Sinh (A : Float) return Float is
------------------------------------------------------------------------------
-- Sinh - returns the hyperbolic sine of its argument
------------------------------------------------------------------------------
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 Float'Machine_Radix = 2 and then Y <= Ln_Float_Large then
Z := Exp (Y);
Z := (Z * 0.5 - 0.5 / Z);
if A < 0.0 then
return -Z;
else
return Z;
end if;
end if;
----Reduce the argument.
if Float'Machine_Radix /= 10 then
W := Y - 8#0.542714#;
else
W := Y - 0.6932;
end if;
if W > Ln_Float_Large then
raise Numeric_Error;
end if;
----Now take the EXP of the argument and compute the result with corrections
-- based on the machine radix.
Z := Exp (W);
if (Float'Machine_Radix = 10 and then
W <= 1.16 * Float (Float'Machine_Mantissa + 1)) then
Z := Z - 0.24997_35916_74870_15965 / Z;
elsif (Float'Machine_Radix /= 10 and then
Float'Machine_Radix /= 2 and then
W <= 0.35 * Float (Float'Machine_Mantissa + 1)) then
Z := Z - 0.24999_30850_04514_99336 / Z;
end if;
if Float'Machine_Radix /= 10 then
Z := Z + 0.13830_27787_96019_02638E-04 * Z;
else
Z := Z + 0.52820_83502_58748_52469E-04 * Z;
end if;
----Return the result with the correct sign.
if A < 0.0 then
return -Z;
else
return Z;
end if;
----Compute the result for "small" arguments here.
else
if Y < Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
return A;
end if;
Y := A * A;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 51) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 51) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 51) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 51) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 16) then
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);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 41) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 41) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 41) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 41) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 13) then
W := (((0.77239_39820_29419_23E-02 * Y +
0.13286_42866_92242_29E+01) * Y +
0.85943_28483_85490_10E+02) * Y +
0.23941_43592_30500_69E+04) * Y
/ ((Y - 0.20258_33686_64278_69E+03) * Y +
0.14364_86155_38302_92E+05);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 8) then
W := ((0.34364_14035_8506E+00 * Y +
0.31359_75645_6058E+02) * Y +
0.10622_28883_7151E+04) * Y
/ ((Y - 0.13051_01250_9199E+03) * Y +
0.63733_73302_1822E+04);
else
W := (-0.19033_3399E+00 * Y - 0.71379_3159E+01) * Y
/ (Y - 0.42827_7109E+02);
end if;
W := A + A * W;
return W;
end if;
end Sinh;
--\f
function Cosh (A : Float) return Float is
------------------------------------------------------------------------------
-- Cosh - returns the hyperbolic cosine of its argument
------------------------------------------------------------------------------
Y : Float := abs A;
W : Float;
Z : Float;
begin
----Base-2 machines can compute things the easy way for "small" arguments.
if (Float'Machine_Radix = 2 and then Y <= Ln_Float_Large) or else
(Float'Machine_Radix /= 2 and then Y <= 1.0) then
Z := Exp (Y);
Z := (Z * 0.5 + 0.5 / Z);
return Z;
end if;
----Reduce the argument.
if Float'Machine_Radix /= 10 then
W := Y - 8#0.542714#;
else
W := Y - 0.6932;
end if;
if W > Ln_Float_Large then
raise Numeric_Error;
end if;
----Now take the EXP of the argument and compute the result with corrections
-- based on the machine radix.
Z := Exp (W);
if (Float'Machine_Radix = 10 and then
W <= 1.16 * Float (Float'Machine_Mantissa + 1)) then
Z := Z + 0.24997_35916_74870_15965 / Z;
elsif (Float'Machine_Radix /= 10 and then
Float'Machine_Radix /= 2 and then
W <= 0.35 * Float (Float'Machine_Mantissa + 1)) then
Z := Z + 0.24999_30850_04514_99336 / Z;
end if;
if Float'Machine_Radix /= 10 then
Z := Z + 0.13830_27787_96019_02638E-04 * Z;
else
Z := Z + 0.52820_83502_58748_52469E-04 * Z;
end if;
----Return the result.
return Z;
end Cosh;
--\f
function Tanh (A : Float) return Float is
------------------------------------------------------------------------------
-- Tanh - returns the hyperbolic tangent of its argument
------------------------------------------------------------------------------
F : Float := abs A;
G : Float;
Result : Float;
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 / (Exp (F + F) + 1.0);
Result := Result + Result;
----Very small arguments result in themselves.
elsif F < Float (Float'Machine_Radix) **
(-Float'Machine_Mantissa / 2) then
Result := F;
----All other arguments get approximated here.
else
G := F * F;
if (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 49) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 49) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 49) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 49) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 15) then
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);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 37) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 37) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 37) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 37) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 12) then
Result := ((-0.36242_42193_46421_73E-03 * G -
0.92318_68945_14261_77E+00) * G -
0.19059_52242_69822_92E+02) * G
/ ((G + 0.25640_98759_51789_75E+02) * G +
0.57178_56728_09658_17E+02);
elsif (Float'Machine_Radix = 2 and then
Float'Machine_Mantissa >= 25) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 >= 25) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 >= 25) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 >= 25) or else
(Float'Machine_Radix = 10 and then
Float'Machine_Mantissa >= 9) then
Result := (-0.93363_47565_2401E+00 * G -
0.21063_95800_0245E+02)
/ ((G + 0.28077_65347_0471E+02) * G +
0.63191_87401_5582E+02);
else
Result := (-0.38310_10665E-02 * G - 0.82377_28127E+00) * G
/ (G + 0.24713_19654E+01);
end if;
Result := F + F * Result;
end if;
----Return a result with the correct sign.
if A < 0.0 then
return -Result;
else
return Result;
end if;
end Tanh;
--\f
function Sech (A : Float) return Float is
begin
return 1.0 / Cosh (A);
end Sech;
--\f
function Csch (A : Float) return Float is
begin
return 1.0 / Sinh (A);
end Csch;
--\f
function Coth (A : Float) return Float is
begin
return 1.0 / Tanh (A);
end Coth;
--\f
function Asinh (A : Float) return Float is
begin
return Ln (A + Sqrt (A * A + 1.0));
end Asinh;
--\f
function Acosh (A : Float) return Float is
begin
return Ln (A + Sqrt (A * A - 1.0));
end Acosh;
--\f
function Atanh (A : Float) return Float is
begin
return Ln ((1.0 + A) / (1.0 - A)) / 2.0;
end Atanh;
--\f
function Asech (A : Float) return Float is
begin
return Acosh (1.0 / A);
end Asech;
--\f
function Acsch (A : Float) return Float is
begin
return Asinh (1.0 / A);
end Acsch;
--\f
function Acoth (A : Float) return Float is
begin
return Atanh (1.0 / A);
end Acoth;
--\f
function Hex (Dig : Integer32) return Character is
------------------------------------------------------------------------------
-- Returns '0'..'9' or 'A'..'F'; Dig in 0..15.
------------------------------------------------------------------------------
begin
if Dig <= 9 then
return Character'Val (Dig + Character'Pos ('0'));
else
return Character'Val (Dig - 10 + Character'Pos ('A'));
end if;
end Hex;
--\f
function Cat (Frac : Float; Bits : Natural32) return String is
------------------------------------------------------------------------------
-- Convert a fraction that is in the range 1.0/Float'Machine_Radix <= X < 1.0
-- into a bunch of Hex digits with an implicit initial decimal point.
-- This routine is only called when Float'Machine_Radix /= 10
------------------------------------------------------------------------------
Valstr : String (1 .. 4);
I : Integer32;
Nxt : Float;
begin
----Get the next 16 bits off of the top of the fraction.
Nxt := Frac * 16.0 ** 4;
I := Integer32 (Nxt);
if Float (I) > Nxt then
I := I - 1;
end if;
----Remove those 16 bits from the fraction. Do it this "strange" way so that
-- we don't get wobbling-precision problems on base-4/8/16 machines.
Nxt := Frac - Float (I) / 16.0 ** 4;
Nxt := Nxt * 16.0 ** 4;
----Convert the 4 hex digits to a string.
Valstr (1) := Hex (I / 16 ** 3);
Valstr (2) := Hex (I / 16 ** 2 rem 16);
Valstr (3) := Hex (I / 16 rem 16);
Valstr (4) := Hex (I rem 16);
----Get the rest of the digits and return.
if Bits > 16 then
return Valstr & Cat (Nxt, Bits - 16);
else
return Valstr;
end if;
end Cat;
--\f
--/ if R1000 then
function Debug_Tools_Float_Image
(Value : Float;
Level : Natural;
Prefix : String;
Expand_Pointers : Boolean) return String is
------------------------------------------------------------------------------
-- Return a string containing the normal base-10 printed representation of the
-- string and (if Machine_Radix /= 10) the Machine_Radix representation.
------------------------------------------------------------------------------
Frac : Float;
Expo : Integer32;
Str : String (1 .. Float'Digits + 12);
Bits : Natural32;
begin
----Get the base 10 string.
Float_Text_Io.Put (Str, Value, Aft => Float'Digits + 3, Exp => 5);
----Compute the number of bits in the mantissa.
case Float'Machine_Radix is
when 2 =>
Bits := Float'Machine_Mantissa;
when 4 =>
Bits := Float'Machine_Mantissa * 2;
when 8 =>
Bits := Float'Machine_Mantissa * 3;
when 16 =>
Bits := Float'Machine_Mantissa * 4;
when others =>
return Str;
end case;
----Get the fraction and the exponent.
Fraction_Exponent (Value, Frac, Expo);
----Create the string containing the Machine_Radix version of the string and
-- return it plus the base-10 string.
return Str & " 0."
& Cat (Frac, Bits)
& " *" & Natural32'Image (Float'Machine_Radix)
& " **" & Integer32'Image (Expo);
end Debug_Tools_Float_Image;
procedure Register_Float is
new Debug_Tools.Register (Float,
Debug_Tools_Float_Image);
--/ end if; -- R1000
--\f
begin
----Register our preferred debug image routine.
--/ if R1000 then
Register_Float;
--/ end if;
----Ln of the largest and of the smallest numbers that exist.
Ln_Float_Large := Ln (Float'Last);
Ln_Float_Small := Ln (2.0 ** (Float'Machine_Emin - 1));
----Used in argument reduciton in Tanh.
Tanh_Xbig :=
(Ln_2 +
Float (Float'Machine_Mantissa + 1) * Ln (Float (Float'Machine_Radix)))
/ 2.0;
----Cot_Eps1 is smallest positive number that can be reciprocated without
-- Numeric_Error.
if Float'Machine_Emax < -Float'Machine_Emin then
Cot_Eps1 := Float (Float'Machine_Radix) ** (-Float'Machine_Emax + 1);
else
Cot_Eps1 := Float (Float'Machine_Radix) ** (Float'Machine_Emin - 1);
end if;
----Set up arrays of values for the Expx routine. We must do this dynamicly
-- because this package is intended to be portable across architectures and
-- Ada provides us with no way to do this staticly. See pages 98 and 106 of
-- the William J. Cody book for details of what and why.
declare
type Natural32_Array is array (Positive range <>) of Natural32;
procedure Expx_Stuff (A1 : in out Float;
A2 : in out Float;
N1 : Natural32_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 : Integer32;
Mask : Integer32;
Shift : Float;
Tmp : Integer;
begin
if Float'Machine_Radix /= 10 then
----Handle all 2**N-base floating points here.
case Float'Machine_Radix is
when 2 =>
Bits := Float'Machine_Mantissa;
when 4 =>
Bits := Float'Machine_Mantissa * 2;
when 8 =>
Bits := Float'Machine_Mantissa * 3;
when 16 =>
Bits := Float'Machine_Mantissa * 4;
when others =>
raise Program_Error;
end case;
if Bits >= 60 then
A1 := Float (N1 (1)) / 2.0 ** 30 +
Float (N1 (2)) / 2.0 ** 60;
Mask := 2 ** Integer (90 - Bits);
A1 := A1 + Float (N1 (3) - N1 (3) rem Mask) / 2.0 ** 90;
A2 := Float (N1 (3) rem Mask) / 2.0 ** 90;
elsif Bits >= 30 then
A1 := Float (N1 (1)) / 2.0 ** 30;
Mask := 2 ** Integer (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;
else
Mask := 2 ** Integer (30 - Bits);
A1 := Float (N1 (1) - N1 (1) rem Mask) / 2.0 ** 30;
A2 := Float (N1 (1) rem Mask) / 2.0 ** 30 +
Float (N1 (2)) / 2.0 ** 60;
end if;
else
----Handle base-10 floating-point here.
Shift := 1.0 / 10.0 ** 5;
Bits := 0;
A1 := 0.0;
for I in 1 .. 6 loop
A1 := A1 + Float (N1 (I)) * Shift;
Shift := Shift / 10.0 ** 5;
Bits := Bits + 5;
if Bits + 5 > Float'Machine_Mantissa then
Tmp := I + 1;
goto Do_A2;
end if;
end loop;
A2 := 0.0;
return;
<<Do_A2>> null;
Bits := Float'Machine_Mantissa - Bits;
Mask := 10 ** Integer (5 - Bits);
A1 := A1 + Float (N1 (Tmp) - N1 (Tmp) rem Mask) * Shift;
A2 := Float (N1 (Tmp) rem Mask) * Shift;
for I in Tmp + 1 .. 6 loop
A2 := A2 + Float (N1 (I)) * Shift;
Shift := Shift / 10.0 ** 5;
end loop;
end if;
end Expx_Stuff;
begin
if Float'Machine_Radix /= 10 then
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#));
else
Expx_A1 (1) := 1.0;
Expx_A1 (3) := 0.63095_73444_80193_24943_43601_366;
Expx_A1 (5) := 0.39810_71705_53497_25077_02523_051;
Expx_A1 (7) := 0.25118_86431_50958_01110_85032_068;
Expx_A1 (9) := 0.15848_93192_46111_34852_02101_373;
Expx_A1 (11) := 0.1;
Expx_Stuff (Expx_A1 (2), Expx_A2 (1),
(79432, 82347, 24281, 50206, 59182, 82800));
Expx_Stuff (Expx_A1 (4), Expx_A2 (2),
(50118, 72336, 27272, 28500, 15541, 86900));
Expx_Stuff (Expx_A1 (6), Expx_A2 (3),
(31622, 77660, 16837, 93319, 98893, 54400));
Expx_Stuff (Expx_A1 (8), Expx_A2 (4),
(19952, 62314, 96887, 96013, 52455, 39700));
Expx_Stuff (Expx_A1 (10), Expx_A2 (5),
(12589, 25411, 79416, 72104, 23954, 10600));
end if;
end;
----If this is not true then this module will be producing values that are
-- accurate to at-most 60-bits/18-digits. There aren't too many machines
-- with larger floating-point than that but we definitely won't function
-- properly (ie. with the desired accuracy). If you don't care about this
-- then comment out the Assert.
if ((Float'Machine_Radix = 2 and then Float'Machine_Mantissa <= 60) or else
(Float'Machine_Radix = 4 and then
Float'Machine_Mantissa * 2 <= 60) or else
(Float'Machine_Radix = 8 and then
Float'Machine_Mantissa * 3 <= 60) or else
(Float'Machine_Radix = 16 and then
Float'Machine_Mantissa * 4 <= 60) or else
(Float'Machine_Radix = 10 and then Float'Machine_Mantissa <= 18)) then
null;
else
raise Program_Error;
end if;
end Trig;