summaryrefslogtreecommitdiff
path: root/libm/ceilfloor.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ceilfloor.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
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
Diffstat (limited to 'libm/ceilfloor.c')
-rw-r--r--libm/ceilfloor.c179
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__ */