/ base / math / ceilfloor.c
ceilfloor.c
  1  #if defined(__ppc__)
  2  /*******************************************************************************
  3  *                                                                              *
  4  *      File ceilfloor.c,                                                       *
  5  *      Function ceil(x) and floor(x),                                          *
  6  *      Implementation of ceil and floor for the PowerPC.                       *
  7  *                                                                              *
  8  *      Copyright � 1991 Apple Computer, Inc.  All rights reserved.             *
  9  *                                                                              *
 10  *      Written by Ali Sazegari, started on November 1991,                      *
 11  *                                                                              *
 12  *      based on math.h, library code for Macintoshes with a 68881/68882        *
 13  *      by Jim Thomas.                                                          *
 14  *                                                                              *
 15  *      W A R N I N G:  This routine expects a 64 bit double model.             *
 16  *                                                                              *
 17  *      December 03 1992: first rs6000 port.                                    *
 18  *      July     14 1993: comment changes and addition of #pragma fenv_access.  *
 19  *	 May      06 1997: port of the ibm/taligent ceil and floor routines.     *
 20  *	 April    11 2001: first port to os x using gcc.				 *
 21  *	 June	    13 2001: replaced __setflm with in-line assembly			 *
 22  *                                                                              *
 23  *******************************************************************************/
 24  
 25  #if !defined(__ppc__)
 26  #define asm(x)
 27  #endif
 28  
 29  static const double        twoTo52  = 4503599627370496.0;
 30  static const unsigned long signMask = 0x80000000ul;
 31  
 32  typedef union
 33        {
 34        struct {
 35  #if defined(__BIG_ENDIAN__)
 36  	unsigned long int hi;
 37  	unsigned long int lo;
 38  #else
 39  	unsigned long int lo;
 40  	unsigned long int hi;
 41  #endif
 42        } words;
 43        double dbl;
 44        } DblInHex;
 45  
 46  /*******************************************************************************
 47  *            Functions needed for the computation.                             *
 48  *******************************************************************************/
 49  
 50  /*******************************************************************************
 51  *      Ceil(x) returns the smallest integer not less than x.                   *
 52  *******************************************************************************/
 53  
 54  double ceil ( double x )
 55  	{
 56  	DblInHex xInHex,OldEnvironment;
 57  	register double y;
 58  	register unsigned long int xhi;
 59  	register int target;
 60  	
 61  	xInHex.dbl = x;
 62  	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
 63  	target = ( xInHex.words.hi < signMask );
 64  	
 65  	if ( xhi < 0x43300000ul ) 
 66  /*******************************************************************************
 67  *      Is |x| < 2.0^52?                                                        *
 68  *******************************************************************************/
 69  		{
 70  		if ( xhi < 0x3ff00000ul ) 
 71  /*******************************************************************************
 72  *      Is |x| < 1.0?                                                           *
 73  *******************************************************************************/
 74  			{
 75  			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
 76  				return ( x );
 77  			else 
 78  				{			                // inexact case
 79  				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 80  				OldEnvironment.words.lo |= 0x02000000ul;
 81  				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 82  				if ( target )
 83  					return ( 1.0 );
 84  				else
 85  					return ( -0.0 );
 86  				}
 87  			}
 88  /*******************************************************************************
 89  *      Is 1.0 < |x| < 2.0^52?                                                  *
 90  *******************************************************************************/
 91  		if ( target ) 
 92  			{
 93  			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
 94  			if ( y < x )
 95  				return ( y + 1.0 );
 96  			else
 97  				return ( y );
 98  			}
 99  		
100  		else 
101  			{
102  			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
103  			if ( y < x )
104  				return ( y + 1.0 );
105  			else
106  				return ( y );
107  			}
108  		}
109  /*******************************************************************************
110  *      |x| >= 2.0^52 or x is a NaN.                                            *
111  *******************************************************************************/
112  	return ( x );
113  	}
114  
115  /*******************************************************************************
116  *      Floor(x) returns the largest integer not greater than x.                *
117  *******************************************************************************/
118  
119  double floor ( double x )
120  	{
121  	DblInHex xInHex,OldEnvironment;
122  	register double y;
123  	register unsigned long int xhi;
124  	register long int target;
125  	
126  	xInHex.dbl = x;
127  	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
128  	target = ( xInHex.words.hi < signMask );
129  	
130  	if ( xhi < 0x43300000ul ) 
131  /*******************************************************************************
132  *      Is |x| < 2.0^52?                                                        *
133  *******************************************************************************/
134  		{
135  		if ( xhi < 0x3ff00000ul ) 
136  /*******************************************************************************
137  *      Is |x| < 1.0?                                                           *
138  *******************************************************************************/
139  			{
140  			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
141  				return ( x );
142  			else 
143  				{			                // inexact case
144  				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
145  				OldEnvironment.words.lo |= 0x02000000ul;
146  				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
147  				if ( target )
148  					return ( 0.0 );
149  				else
150  					return ( -1.0 );
151  				}
152  			}
153  /*******************************************************************************
154  *      Is 1.0 < |x| < 2.0^52?                                                  *
155  *******************************************************************************/
156  		if ( target ) 
157  			{
158  			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
159  			if ( y > x )
160  				return ( y - 1.0 );
161  			else
162  				return ( y );
163  			}
164  		
165  		else 
166  			{
167  			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
168  			if ( y > x )
169  				return ( y - 1.0 );
170  			else
171  				return ( y );
172  			}
173  		}
174  /*******************************************************************************
175  *      |x| >= 2.0^52 or x is a NaN.                                            *
176  *******************************************************************************/
177  	return ( x );
178  	}
179  #endif /* __ppc__ */