split out nearbyint, round, trunc from libm/powerpc/s_modf.c

5 files changed, 238 insertions, 197 deletions

diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile
index 9e4a309da..342c182ba 100644
--- a/libm/powerpc/Makefile
+++ b/libm/powerpc/Makefile
@@ -38,7 +38,8 @@ LIBM=../libm.a
 CFLAGS+=-D_IEEE_LIBM -D_ISOC99_SOURCE -D_SVID_SOURCE
 
 ifeq ($(strip $(DO_C99_MATH)),y)
-CSRC = s_ceil.c s_floor.c s_ldexp.c s_frexp.c s_logb.c s_modf.c w_scalb.c s_copysign.c s_rint.c
+CSRC = s_ceil.c s_copysign.c s_floor.c s_frexp.c s_ldexp.c s_logb.c s_modf.c \
+	s_nearbyint.c s_rint.c s_round.c s_trunc.c w_scalb.c
 else
 CSRC =
 endif
diff --git a/libm/powerpc/s_modf.c b/libm/powerpc/s_modf.c
index 403c54b79..f4344bda8 100644
--- a/libm/powerpc/s_modf.c
+++ b/libm/powerpc/s_modf.c
@@ -4,9 +4,8 @@
 **      Contains: C source code for implementations of floating-point
 **                functions which round to integral value or format, as
 **                defined in header <fp.h>.  In particular, this file
-**                contains implementations of functions rint, nearbyint,
-**                rinttol, round, roundtol, trunc, modf and modfl.  This file
-**                targets PowerPC or Power platforms.
+**                contains implementations of functions rinttol, roundtol,
+**                modf and modfl.  This file targets PowrPC or Power platforms.
 **
 **      Written by: A. Sazegari, Apple AltiVec Group
 **	    Created originally by Jon Okada, Apple Numerics Group
@@ -66,44 +65,11 @@ typedef union
 static const unsigned long int signMask = 0x80000000ul;
 static const double twoTo52      = 4503599627370496.0;
 static const double doubleToLong = 4503603922337792.0;	            // 2^52
-static const DblInHex Huge       = {{ 0x7FF00000, 0x00000000 }};
 static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
 
 
 /*******************************************************************************
 *                                                                              *
-*     The function nearbyint rounds its double argument to integral value      *
-*     according to the current rounding direction and returns the result in    *
-*     double format.  This function does not signal inexact.                   *
-*                                                                              *
-********************************************************************************
-*                                                                              *
-*     This function calls fabs and copysign.		                         *
-*                                                                              *
-*******************************************************************************/
-
-double nearbyint ( double x )
-      {
-	double y;
-	double OldEnvironment;
-
-	y = twoTo52;
-
-	asm ("mffs %0" : "=f" (OldEnvironment));	/* get the environement */
-
-      if ( fabs ( x ) >= y )                          /* huge case is exact */
-            return x;
-      if ( x < 0 ) y = -y;                                   /* negative case */
-      y = ( x + y ) - y;                                    /* force rounding */
-      if ( y == 0.0 )                        /* zero results mirror sign of x */
-            y = copysign ( y, x );
-//	restore old flags
-	asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
-      return ( y );
-	}
-
-/*******************************************************************************
-*                                                                              *
 *     The function rinttol converts its double argument to integral value      *
 *     according to the current rounding direction and returns the result in    *
 *     long int format.  This conversion signals invalid if the argument is a   *
@@ -197,99 +163,6 @@ long int rinttol ( double x )
 
 /*******************************************************************************
 *                                                                              *
-*     The function round rounds its double argument to integral value          *
-*     according to the "add half to the magnitude and truncate" rounding of    *
-*     Pascal's Round function and FORTRAN's ANINT function and returns the     *
-*     result in double format.  This function signals inexact if an ordered    *
-*     return value is not equal to the operand.                                *
-*                                                                              *
-*******************************************************************************/
-
-double round ( double x )
-      {
-      DblInHex argument, OldEnvironment;
-      register double y, z;
-      register unsigned long int xHead;
-      register long int target;
-
-      argument.dbl = x;
-      xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
-      target = ( argument.words.hi < signMask );     // flag positive sign
-
-      if ( xHead < 0x43300000ul )
-/*******************************************************************************
-*     Is |x| < 2.0^52?                                                        *
-*******************************************************************************/
-            {
-            if ( xHead < 0x3ff00000ul )
-/*******************************************************************************
-*     Is |x| < 1.0?                                                           *
-*******************************************************************************/
-                  {
-			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
-                  if ( xHead < 0x3fe00000ul )
-/*******************************************************************************
-*     Is |x| < 0.5?                                                           *
-*******************************************************************************/
-                        {
-                        if ( ( xHead | argument.words.lo ) != 0ul )
-                              OldEnvironment.words.lo |= 0x02000000ul;
-				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
-                        if ( target )
-                              return ( 0.0 );
-                        else
-                              return ( -0.0 );
-                        }
-/*******************************************************************************
-*     Is 0.5 ² |x| < 1.0?                                                      *
-*******************************************************************************/
-                  OldEnvironment.words.lo |= 0x02000000ul;
-			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
-                  if ( target )
-                        return ( 1.0 );
-                  else
-                        return ( -1.0 );
-                  }
-/*******************************************************************************
-*     Is 1.0 < |x| < 2.0^52?                                                   *
-*******************************************************************************/
-            if ( target )
-                  {                                     // positive x
-                  y = ( x + twoTo52 ) - twoTo52;        // round at binary point
-                  if ( y == x )                         // exact case
-                        return ( x );
-                  z = x + 0.5;                          // inexact case
-                  y = ( z + twoTo52 ) - twoTo52;        // round at binary point
-                  if ( y > z )
-                        return ( y - 1.0 );
-                  else
-                        return ( y );
-                  }
-
-/*******************************************************************************
-*     Is x < 0?                                                                *
-*******************************************************************************/
-            else
-                  {
-                  y = ( x - twoTo52 ) + twoTo52;        // round at binary point
-                  if ( y == x )
-                        return ( x );
-                  z = x - 0.5;
-                  y = ( z - twoTo52 ) + twoTo52;        // round at binary point
-                  if ( y < z )
-                        return ( y + 1.0 );
-                  else
-                  return ( y );
-                  }
-            }
-/*******************************************************************************
-*      |x| >= 2.0^52 or x is a NaN.                                            *
-*******************************************************************************/
-      return ( x );
-      }
-
-/*******************************************************************************
-*                                                                              *
 *     The function roundtol converts its double argument to integral format    *
 *     according to the "add half to the magnitude and chop" rounding mode of   *
 *     Pascal's Round function and FORTRAN's NINT function.  This conversion    *
@@ -392,73 +265,6 @@ long int roundtol ( double x )
 	}
 
 /*******************************************************************************
-*                                                                              *
-*     The function trunc truncates its double argument to integral value       *
-*     and returns the result in double format.  This function signals          *
-*     inexact if an ordered return value is not equal to the operand.          *
-*                                                                              *
-*******************************************************************************/
-
-double trunc ( double x )
-      {
-	DblInHex argument,OldEnvironment;
-	register double y;
-	register unsigned long int xhi;
-	register long int target;
-
-	argument.dbl = x;
-	xhi = argument.words.hi & 0x7fffffffUL;	      	// xhi <- high half of |x|
-	target = ( argument.words.hi < signMask );	      	// flag positive sign
-
-	if ( xhi < 0x43300000ul )
-/*******************************************************************************
-*     Is |x| < 2.0^53?                                                         *
-*******************************************************************************/
-		{
-		if ( xhi < 0x3ff00000ul )
-/*******************************************************************************
-*     Is |x| < 1.0?                                                            *
-*******************************************************************************/
-			{
-			if ( ( xhi | argument.words.lo ) != 0ul )
-				{                             	// raise deserved INEXACT
-				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
-				OldEnvironment.words.lo |= 0x02000000ul;
-				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
-				}
-			if ( target )	                  	// return properly signed zero
-				return ( 0.0 );
-			else
-				return ( -0.0 );
-			}
-/*******************************************************************************
-*     Is 1.0 < |x| < 2.0^52?                                                   *
-*******************************************************************************/
-		if ( target )
-			{
-			y = ( x + twoTo52 ) - twoTo52;      	// round at binary point
-			if ( y > x )
-				return ( y - 1.0 );
-			else
-				return ( y );
-			}
-
-		else
-			{
-			y = ( x - twoTo52 ) + twoTo52;      	// round at binary point.
-			if ( y < x )
-				return ( y + 1.0 );
-			else
-				return ( y );
-			}
-		}
-/*******************************************************************************
-*      Is |x| >= 2.0^52 or x is a NaN.                                         *
-*******************************************************************************/
-	return ( x );
-	}
-
-/*******************************************************************************
 *     The modf family of functions separate a floating-point number into its   *
 *     fractional and integral parts, returning the fractional part and writing *
 *     the integral part in floating-point format to the object pointed to by a *
diff --git a/libm/powerpc/s_nearbyint.c b/libm/powerpc/s_nearbyint.c
new file mode 100644
index 000000000..f2d7ded35
--- /dev/null
+++ b/libm/powerpc/s_nearbyint.c
@@ -0,0 +1,36 @@
+#include <limits.h>
+#include <math.h>
+
+/*******************************************************************************
+*                                                                              *
+*     The function nearbyint rounds its double argument to integral value      *
+*     according to the current rounding direction and returns the result in    *
+*     double format.  This function does not signal inexact.                   *
+*                                                                              *
+********************************************************************************
+*                                                                              *
+*     This function calls fabs and copysign.		                         *
+*                                                                              *
+*******************************************************************************/
+
+static const double twoTo52      = 4503599627370496.0;
+
+double nearbyint ( double x )
+      {
+	double y;
+	double OldEnvironment;
+
+	y = twoTo52;
+
+	asm ("mffs %0" : "=f" (OldEnvironment));	/* get the environement */
+
+      if ( fabs ( x ) >= y )                          /* huge case is exact */
+            return x;
+      if ( x < 0 ) y = -y;                                   /* negative case */
+      y = ( x + y ) - y;                                    /* force rounding */
+      if ( y == 0.0 )                        /* zero results mirror sign of x */
+            y = copysign ( y, x );
+//	restore old flags
+	asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
+      return ( y );
+	}
diff --git a/libm/powerpc/s_round.c b/libm/powerpc/s_round.c
new file mode 100644
index 000000000..81f4d0fec
--- /dev/null
+++ b/libm/powerpc/s_round.c
@@ -0,0 +1,112 @@
+#include <limits.h>
+#include <math.h>
+
+typedef union
+      {
+      struct {
+#if defined(__BIG_ENDIAN__)
+        unsigned long int hi;
+        unsigned long int lo;
+#else
+        unsigned long int lo;
+        unsigned long int hi;
+#endif
+      } words;
+      double dbl;
+      } DblInHex;
+
+static const unsigned long int signMask = 0x80000000ul;
+static const double twoTo52      = 4503599627370496.0;
+
+/*******************************************************************************
+*                                                                              *
+*     The function round rounds its double argument to integral value          *
+*     according to the "add half to the magnitude and truncate" rounding of    *
+*     Pascal's Round function and FORTRAN's ANINT function and returns the     *
+*     result in double format.  This function signals inexact if an ordered    *
+*     return value is not equal to the operand.                                *
+*                                                                              *
+*******************************************************************************/
+
+double round ( double x )
+      {
+      DblInHex argument, OldEnvironment;
+      register double y, z;
+      register unsigned long int xHead;
+      register long int target;
+
+      argument.dbl = x;
+      xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
+      target = ( argument.words.hi < signMask );     // flag positive sign
+
+      if ( xHead < 0x43300000ul )
+/*******************************************************************************
+*     Is |x| < 2.0^52?                                                        *
+*******************************************************************************/
+            {
+            if ( xHead < 0x3ff00000ul )
+/*******************************************************************************
+*     Is |x| < 1.0?                                                           *
+*******************************************************************************/
+                  {
+			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
+                  if ( xHead < 0x3fe00000ul )
+/*******************************************************************************
+*     Is |x| < 0.5?                                                           *
+*******************************************************************************/
+                        {
+                        if ( ( xHead | argument.words.lo ) != 0ul )
+                              OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+                        if ( target )
+                              return ( 0.0 );
+                        else
+                              return ( -0.0 );
+                        }
+/*******************************************************************************
+*     Is 0.5 ² |x| < 1.0?                                                      *
+*******************************************************************************/
+                  OldEnvironment.words.lo |= 0x02000000ul;
+			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+                  if ( target )
+                        return ( 1.0 );
+                  else
+                        return ( -1.0 );
+                  }
+/*******************************************************************************
+*     Is 1.0 < |x| < 2.0^52?                                                   *
+*******************************************************************************/
+            if ( target )
+                  {                                     // positive x
+                  y = ( x + twoTo52 ) - twoTo52;        // round at binary point
+                  if ( y == x )                         // exact case
+                        return ( x );
+                  z = x + 0.5;                          // inexact case
+                  y = ( z + twoTo52 ) - twoTo52;        // round at binary point
+                  if ( y > z )
+                        return ( y - 1.0 );
+                  else
+                        return ( y );
+                  }
+
+/*******************************************************************************
+*     Is x < 0?                                                                *
+*******************************************************************************/
+            else
+                  {
+                  y = ( x - twoTo52 ) + twoTo52;        // round at binary point
+                  if ( y == x )
+                        return ( x );
+                  z = x - 0.5;
+                  y = ( z - twoTo52 ) + twoTo52;        // round at binary point
+                  if ( y < z )
+                        return ( y + 1.0 );
+                  else
+                  return ( y );
+                  }
+            }
+/*******************************************************************************
+*      |x| >= 2.0^52 or x is a NaN.                                            *
+*******************************************************************************/
+      return ( x );
+      }
diff --git a/libm/powerpc/s_trunc.c b/libm/powerpc/s_trunc.c
new file mode 100644
index 000000000..4b61355ea
--- /dev/null
+++ b/libm/powerpc/s_trunc.c
@@ -0,0 +1,86 @@
+#include <limits.h>
+#include <math.h>
+
+typedef union
+      {
+      struct {
+#if defined(__BIG_ENDIAN__)
+        unsigned long int hi;
+        unsigned long int lo;
+#else
+        unsigned long int lo;
+        unsigned long int hi;
+#endif
+      } words;
+      double dbl;
+      } DblInHex;
+
+static const unsigned long int signMask = 0x80000000ul;
+static const double twoTo52      = 4503599627370496.0;
+
+/*******************************************************************************
+*                                                                              *
+*     The function trunc truncates its double argument to integral value       *
+*     and returns the result in double format.  This function signals          *
+*     inexact if an ordered return value is not equal to the operand.          *
+*                                                                              *
+*******************************************************************************/
+
+double trunc ( double x )
+      {
+	DblInHex argument,OldEnvironment;
+	register double y;
+	register unsigned long int xhi;
+	register long int target;
+
+	argument.dbl = x;
+	xhi = argument.words.hi & 0x7fffffffUL;	      	// xhi <- high half of |x|
+	target = ( argument.words.hi < signMask );	      	// flag positive sign
+
+	if ( xhi < 0x43300000ul )
+/*******************************************************************************
+*     Is |x| < 2.0^53?                                                         *
+*******************************************************************************/
+		{
+		if ( xhi < 0x3ff00000ul )
+/*******************************************************************************
+*     Is |x| < 1.0?                                                            *
+*******************************************************************************/
+			{
+			if ( ( xhi | argument.words.lo ) != 0ul )
+				{                             	// raise deserved INEXACT
+				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
+				OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+				}
+			if ( target )	                  	// return properly signed zero
+				return ( 0.0 );
+			else
+				return ( -0.0 );
+			}
+/*******************************************************************************
+*     Is 1.0 < |x| < 2.0^52?                                                   *
+*******************************************************************************/
+		if ( target )
+			{
+			y = ( x + twoTo52 ) - twoTo52;      	// round at binary point
+			if ( y > x )
+				return ( y - 1.0 );
+			else
+				return ( y );
+			}
+
+		else
+			{
+			y = ( x - twoTo52 ) + twoTo52;      	// round at binary point.
+			if ( y < x )
+				return ( y + 1.0 );
+			else
+				return ( y );
+			}
+		}
+/*******************************************************************************
+*      Is |x| >= 2.0^52 or x is a NaN.                                         *
+*******************************************************************************/
+	return ( x );
+	}