Index: a-numaux-vxworks.ads
===================================================================
--- a-numaux-vxworks.ads	(revision 247135)
+++ a-numaux-vxworks.ads	(working copy)
@@ -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
Index: a-numaux-x86.adb
===================================================================
--- a-numaux-x86.adb	(revision 247135)
+++ a-numaux-x86.adb	(working copy)
@@ -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;
 
    ----------
Index: a-numaux-x86.ads
===================================================================
--- a-numaux-x86.ads	(revision 247135)
+++ a-numaux-x86.ads	(working copy)
@@ -30,7 +30,8 @@
 --                                                                          --
 ------------------------------------------------------------------------------
 
---  Version for the x86, using 64-bit IEEE format with inline asm statements
+--  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;
Index: a-numaux-darwin.adb
===================================================================
--- a-numaux-darwin.adb	(revision 247135)
+++ a-numaux-darwin.adb	(working copy)
@@ -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;
 
    ---------
Index: a-numaux-darwin.ads
===================================================================
--- a-numaux-darwin.ads	(revision 247135)
+++ a-numaux-darwin.ads	(working copy)
@@ -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
Index: a-numaux.ads
===================================================================
--- a-numaux.ads	(revision 247135)
+++ a-numaux.ads	(working copy)
@@ -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 OS/X (special handling of sin/cos for accuracy)
+--    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 64-bit IEEE (inline asm statements)
+--    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
Index: a-numaux-libc-x86.ads
===================================================================
--- a-numaux-libc-x86.ads	(revision 247135)
+++ a-numaux-libc-x86.ads	(working copy)
@@ -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.
