libm/powerpc/classic/s_ceil.c

   1 /*******************************************************************************
   2 *                                                                              *
   3 *      File ceilfloor.c,                                                       *
   4 *      Function ceil(x) and floor(x),                                          *
   5 *      Implementation of ceil and floor for the PowerPC.                       *
   6 *                                                                              *
   7 *      Copyright © 1991 Apple Computer, Inc.  All rights reserved.             *
   8 *                                                                              *
   9 *      Written by Ali Sazegari, started on November 1991,                      *
  10 *                                                                              *
  11 *      based on math.h, library code for Macintoshes with a 68881/68882        *
  12 *      by Jim Thomas.                                                          *
  13 *                                                                              *
  14 *      W A R N I N G:  This routine expects a 64 bit double model.             *
  15 *                                                                              *
  16 *      December 03 1992: first rs6000 port.                                    *
  17 *      July     14 1993: comment changes and addition of #pragma fenv_access.  *
  18 *        May      06 1997: port of the ibm/taligent ceil and floor routines.     *
  19 *        April    11 2001: first port to os x using gcc.                                 *
  20 *        June       13 2001: replaced __setflm with in-line assembly                     *
  21 *                                                                              *
  22 *******************************************************************************/
  23
  24 #include <math.h>
  25 #include <endian.h>
  26
  27 static const double        twoTo52  = 4503599627370496.0;
  28 static const unsigned long signMask = 0x80000000ul;
  29
  30 typedef union
  31       {
  32       struct {
  33 #if (__BYTE_ORDER == __BIG_ENDIAN)
  34         unsigned long int hi;
  35         unsigned long int lo;
  36 #else
  37         unsigned long int lo;
  38         unsigned long int hi;
  39 #endif
  40       } words;
  41       double dbl;
  42       } DblInHex;
  43
  44 /*******************************************************************************
  45 *            Functions needed for the computation.                             *
  46 *******************************************************************************/
  47
  48 /*******************************************************************************
  49 *      Ceil(x) returns the smallest integer not less than x.                   *
  50 *******************************************************************************/
  51
  52 libm_hidden_proto(ceil)
  53 double ceil ( double x )
  54         {
  55         DblInHex xInHex,OldEnvironment;
  56         register double y;
  57         register unsigned long int xhi;
  58         register int target;
  59
  60         xInHex.dbl = x;
  61         xhi = xInHex.words.hi & 0x7fffffffUL;     // xhi is the high half of |x|
  62         target = ( xInHex.words.hi < signMask );
  63
  64         if ( xhi < 0x43300000ul )
  65 /*******************************************************************************
  66 *      Is |x| < 2.0^52?                                                        *
  67 *******************************************************************************/
  68                 {
  69                 if ( xhi < 0x3ff00000ul )
  70 /*******************************************************************************
  71 *      Is |x| < 1.0?                                                           *
  72 *******************************************************************************/
  73                         {
  74                         if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
  75                                 return ( x );
  76                         else
  77                                 {                                       // inexact case
  78                                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));
  79                                 OldEnvironment.words.lo |= 0x02000000ul;
  80                                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  81                                 if ( target )
  82                                         return ( 1.0 );
  83                                 else
  84                                         return ( -0.0 );
  85                                 }
  86                         }
  87 /*******************************************************************************
  88 *      Is 1.0 < |x| < 2.0^52?                                                  *
  89 *******************************************************************************/
  90                 if ( target )
  91                         {
  92                         y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
  93                         if ( y < x )
  94                                 return ( y + 1.0 );
  95                         else
  96                                 return ( y );
  97                         }
  98
  99                 else
 100                         {
 101                         y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
 102                         if ( y < x )
 103                                 return ( y + 1.0 );
 104                         else
 105                                 return ( y );
 106                         }
 107                 }
 108 /*******************************************************************************
 109 *      |x| >= 2.0^52 or x is a NaN.                                            *
 110 *******************************************************************************/
 111         return ( x );
 112         }
 113 libm_hidden_def(ceil)