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

1 files changed, 112 insertions, 0 deletions

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 );
+      }