===================================================================
@@ -36,7 +36,7 @@
package Ada.Numerics.Aux is
pragma Pure;
- type Double is digits 15;
+ type Double is new Long_Float;
-- Type Double is the type used to call the C routines
-- We import these functions directly from C. Note that we label them
===================================================================
@@ -49,8 +49,11 @@
-- for values of Y in the open interval (-0.25, 0.25)
procedure Reduce (X : in out Double; Q : out Natural);
- -- Implements reduction of X by Pi/2. Q is the quadrant of the final
- -- result in the range 0 .. 3. The absolute value of X is at most Pi.
+ -- Implement reduction of X by Pi/2. Q is the quadrant of the final
+ -- result in the range 0..3. The absolute value of X is at most Pi/4.
+ -- It is needed to avoid a loss of accuracy for sin near Pi and cos
+ -- near Pi/2 due to the use of an insufficiently precise value of Pi
+ -- in the range reduction.
pragma Inline (Is_Nan);
pragma Inline (Reduce);
@@ -117,7 +120,7 @@
begin
-- The IEEE NaN values are the only ones that do not equal themselves
- return not (X = X);
+ return X /= X;
end Is_Nan;
---------
@@ -154,32 +157,36 @@
P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
- P4, HM);
P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
- K : Double := X * Two_Over_Pi;
+ K : Double;
+ R : Integer;
+
begin
- -- For X < 2.0**32, all products below are computed exactly.
+ -- For X < 2.0**HM, all products below are computed exactly.
-- Due to cancellation effects all subtractions are exact as well.
-- As no double extended floating-point number has more than 75
-- zeros after the binary point, the result will be the correctly
-- rounded result of X - K * (Pi / 2.0).
+ K := X * Two_Over_Pi;
while abs K >= 2.0**HM loop
K := K * M - (K * M - K);
- X := (((((X - K * P1) - K * P2) - K * P3)
- - K * P4) - K * P5) - K * P6;
+ X :=
+ (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
K := X * Two_Over_Pi;
end loop;
- if K /= K then
+ -- If K is not a number (because X was not finite) raise exception
- -- K is not a number, because X was not finite
-
+ if Is_Nan (K) then
raise Constraint_Error;
end if;
- K := Double'Rounding (K);
- Q := Integer (K) mod 4;
- X := (((((X - K * P1) - K * P2) - K * P3)
- - K * P4) - K * P5) - K * P6;
+ -- Go through an integer temporary so as to use machine instructions
+
+ R := Integer (Double'Rounding (K));
+ Q := R mod 4;
+ K := Double (R);
+ X := (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
end Reduce;
----------
===================================================================
@@ -30,7 +30,8 @@
-- --
------------------------------------------------------------------------------
+-- This version is for the x86 using the 80-bit x86 long double format with
+-- inline asm statements.
package Ada.Numerics.Aux is
pragma Pure;
===================================================================
@@ -36,11 +36,17 @@
-- Local subprograms --
-----------------------
+ function Is_Nan (X : Double) return Boolean;
+ -- Return True iff X is a IEEE NaN value
+
procedure Reduce (X : in out Double; Q : out Natural);
- -- Implements reduction of X by Pi/2. Q is the quadrant of the final
- -- result in the range 0 .. 3. The absolute value of X is at most Pi/4.
+ -- Implement reduction of X by Pi/2. Q is the quadrant of the final
+ -- result in the range 0..3. The absolute value of X is at most Pi/4.
+ -- It is needed to avoid a loss of accuracy for sin near Pi and cos
+ -- near Pi/2 due to the use of an insufficiently precise value of Pi
+ -- in the range reduction.
- -- The following three functions implement Chebishev approximations
+ -- The following two functions implement Chebishev approximations
-- of the trigonometric functions in their reduced domain.
-- These approximations have been computed using Maple.
@@ -51,6 +57,10 @@
pragma Inline (Sine_Approx);
pragma Inline (Cosine_Approx);
+ -------------------
+ -- Cosine_Approx --
+ -------------------
+
function Cosine_Approx (X : Double) return Double is
XX : constant Double := X * X;
begin
@@ -63,6 +73,10 @@
- 16#3.655E64869ECCE#E-14 + 1.0;
end Cosine_Approx;
+ -----------------
+ -- Sine_Approx --
+ -----------------
+
function Sine_Approx (X : Double) return Double is
XX : constant Double := X * X;
begin
@@ -75,6 +89,17 @@
end Sine_Approx;
------------
+ -- Is_Nan --
+ ------------
+
+ function Is_Nan (X : Double) return Boolean is
+ begin
+ -- The IEEE NaN values are the only ones that do not equal themselves
+
+ return X /= X;
+ end Is_Nan;
+
+ ------------
-- Reduce --
------------
@@ -92,6 +117,7 @@
- P4, HM);
P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
K : Double;
+ R : Integer;
begin
-- For X < 2.0**HM, all products below are computed exactly.
@@ -101,7 +127,7 @@
-- rounded result of X - K * (Pi / 2.0).
K := X * Two_Over_Pi;
- while abs K >= 2.0 ** HM loop
+ while abs K >= 2.0**HM loop
K := K * M - (K * M - K);
X :=
(((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
@@ -110,14 +136,16 @@
-- If K is not a number (because X was not finite) raise exception
- if K /= K then
+ if Is_Nan (K) then
raise Constraint_Error;
end if;
- K := Double'Rounding (K);
- Q := Integer (K) mod 4;
- X := (((((X - K * P1) - K * P2) - K * P3)
- - K * P4) - K * P5) - K * P6;
+ -- Go through an integer temporary so as to use machine instructions
+
+ R := Integer (Double'Rounding (K));
+ Q := R mod 4;
+ K := Double (R);
+ X := (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
end Reduce;
---------
===================================================================
@@ -39,7 +39,7 @@
pragma Linker_Options ("-lm");
- type Double is digits 15;
+ type Double is new Long_Float;
-- Type Double is the type used to call the C routines
-- The following functions have been implemented in Ada, since
===================================================================
@@ -40,9 +40,10 @@
-- This version here is for use with normal Unix math functions. Alternative
-- versions are provided for special situations:
+-- a-numaux-darwin For PowerPC/Darwin (special handling of sin/cos)
-- a-numaux-libc-x86 For the x86, using 80-bit long double format
+-- a-numaux-x86 For the x86, using 80-bit long double format with
+-- inline asm statements
-- a-numaux-vxworks For use on VxWorks (where we have no libm.a library)
package Ada.Numerics.Aux is
@@ -50,7 +51,7 @@
pragma Linker_Options ("-lm");
- type Double is digits 15;
+ type Double is new Long_Float;
-- Type Double is the type used to call the C routines
-- We import these functions directly from C. Note that we label them
===================================================================
@@ -37,7 +37,7 @@
pragma Linker_Options ("-lm");
- type Double is digits 18;
+ type Double is new Long_Long_Float;
-- We import these functions directly from C. Note that we label them
-- all as pure functions, because indeed all of them are in fact pure.