Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 179 insertions, 0 deletions

diff --git a/libm/ceilfloor.c b/libm/ceilfloor.c
new file mode 100644
index 000000000..9607435c3
--- /dev/null
+++ b/libm/ceilfloor.c
@@ -0,0 +1,179 @@
+#if defined(__ppc__)
+/*******************************************************************************
+*                                                                              *
+*      File ceilfloor.c,                                                       *
+*      Function ceil(x) and floor(x),                                          *
+*      Implementation of ceil and floor for the PowerPC.                       *
+*                                                                              *
+*      Copyright © 1991 Apple Computer, Inc.  All rights reserved.             *
+*                                                                              *
+*      Written by Ali Sazegari, started on November 1991,                      *
+*                                                                              *
+*      based on math.h, library code for Macintoshes with a 68881/68882        *
+*      by Jim Thomas.                                                          *
+*                                                                              *
+*      W A R N I N G:  This routine expects a 64 bit double model.             *
+*                                                                              *
+*      December 03 1992: first rs6000 port.                                    *
+*      July     14 1993: comment changes and addition of #pragma fenv_access.  *
+*	 May      06 1997: port of the ibm/taligent ceil and floor routines.     *
+*	 April    11 2001: first port to os x using gcc.				 *
+*	 June	    13 2001: replaced __setflm with in-line assembly			 *
+*                                                                              *
+*******************************************************************************/
+
+#if !defined(__ppc__)
+#define asm(x)
+#endif
+
+static const double        twoTo52  = 4503599627370496.0;
+static const unsigned long signMask = 0x80000000ul;
+
+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;
+
+/*******************************************************************************
+*            Functions needed for the computation.                             *
+*******************************************************************************/
+
+/*******************************************************************************
+*      Ceil(x) returns the smallest integer not less than x.                   *
+*******************************************************************************/
+
+double ceil ( double x )
+	{
+	DblInHex xInHex,OldEnvironment;
+	register double y;
+	register unsigned long int xhi;
+	register int target;
+	
+	xInHex.dbl = x;
+	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
+	target = ( xInHex.words.hi < signMask );
+	
+	if ( xhi < 0x43300000ul ) 
+/*******************************************************************************
+*      Is |x| < 2.0^52?                                                        *
+*******************************************************************************/
+		{
+		if ( xhi < 0x3ff00000ul ) 
+/*******************************************************************************
+*      Is |x| < 1.0?                                                           *
+*******************************************************************************/
+			{
+			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
+				return ( x );
+			else 
+				{			                // inexact case
+				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
+				OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+				if ( target )
+					return ( 1.0 );
+				else
+					return ( -0.0 );
+				}
+			}
+/*******************************************************************************
+*      Is 1.0 < |x| < 2.0^52?                                                  *
+*******************************************************************************/
+		if ( target ) 
+			{
+			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
+			if ( y < x )
+				return ( y + 1.0 );
+			else
+				return ( y );
+			}
+		
+		else 
+			{
+			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
+			if ( y < x )
+				return ( y + 1.0 );
+			else
+				return ( y );
+			}
+		}
+/*******************************************************************************
+*      |x| >= 2.0^52 or x is a NaN.                                            *
+*******************************************************************************/
+	return ( x );
+	}
+
+/*******************************************************************************
+*      Floor(x) returns the largest integer not greater than x.                *
+*******************************************************************************/
+
+double floor ( double x )
+	{
+	DblInHex xInHex,OldEnvironment;
+	register double y;
+	register unsigned long int xhi;
+	register long int target;
+	
+	xInHex.dbl = x;
+	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
+	target = ( xInHex.words.hi < signMask );
+	
+	if ( xhi < 0x43300000ul ) 
+/*******************************************************************************
+*      Is |x| < 2.0^52?                                                        *
+*******************************************************************************/
+		{
+		if ( xhi < 0x3ff00000ul ) 
+/*******************************************************************************
+*      Is |x| < 1.0?                                                           *
+*******************************************************************************/
+			{
+			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
+				return ( x );
+			else 
+				{			                // inexact case
+				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
+				OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+				if ( target )
+					return ( 0.0 );
+				else
+					return ( -1.0 );
+				}
+			}
+/*******************************************************************************
+*      Is 1.0 < |x| < 2.0^52?                                                  *
+*******************************************************************************/
+		if ( target ) 
+			{
+			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
+			if ( y > x )
+				return ( y - 1.0 );
+			else
+				return ( y );
+			}
+		
+		else 
+			{
+			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
+			if ( y > x )
+				return ( y - 1.0 );
+			else
+				return ( y );
+			}
+		}
+/*******************************************************************************
+*      |x| >= 2.0^52 or x is a NaN.                                            *
+*******************************************************************************/
+	return ( x );
+	}
+#endif /* __ppc__ */