gonzui

    1: /* @(#)s_ceil.c 5.1 93/09/24 */
    2: /*
    3:  * ====================================================
    4:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    5:  *
    6:  * Developed at SunPro, a Sun Microsystems, Inc. business.
    7:  * Permission to use, copy, modify, and distribute this
    8:  * software is freely granted, provided that this notice 
    9:  * is preserved.
   10:  * ====================================================
   11:  */
   12: 
   13: /* LINTLIBRARY */
   14: 
   15: /*
   16:  * ceil(x)
   17:  * Return x rounded toward -inf to integral value
   18:  * Method:
   19:  *      Bit twiddling.
   20:  * Exception:
   21:  *      Inexact flag raised if x not equal to ceil(x).
   22:  */
   23: 
   24: #include <sys/cdefs.h>
   25: #include <float.h>
   26: #include <math.h>
   27: 
   28: #include "math_private.h"
   29: 
   30: static const double huge = 1.0e300;
   31: 
   32: double
   33: ceil(double x)
   34: {
   35:         int32_t i0,i1,jj0;
   36:         u_int32_t i,j;
   37:         EXTRACT_WORDS(i0,i1,x);
   38:         jj0 = ((i0>>20)&0x7ff)-0x3ff;
   39:         if(jj0<20) {
   40:             if(jj0<0) {        /* raise inexact if x != 0 */
   41:                 if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
   42:                     if(i0<0) {i0=0x80000000;i1=0;} 
   43:                     else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
   44:                 }
   45:             } else {
   46:                 i = (0x000fffff)>>jj0;
   47:                 if(((i0&i)|i1)==0) return x; /* x is integral */
   48:                 if(huge+x>0.0) {      /* raise inexact flag */
   49:                     if(i0>0) i0 += (0x00100000)>>jj0;
   50:                     i0 &= (~i); i1=0;
   51:                 }
   52:             }
   53:         } else if (jj0>51) {
   54:             if(jj0==0x400) return x+x; /* inf or NaN */
   55:             else return x;             /* x is integral */
   56:         } else {
   57:             i = ((u_int32_t)(0xffffffff))>>(jj0-20);
   58:             if((i1&i)==0) return x;    /* x is integral */
   59:             if(huge+x>0.0) {           /* raise inexact flag */
   60:                 if(i0>0) {
   61:                     if(jj0==20) i0+=1; 
   62:                     else {
   63:                         j = i1 + (1<<(52-jj0));
   64:                         if(j<i1) i0+=1;      /* got a carry */
   65:                         i1 = j;
   66:                     }
   67:                 }
   68:                 i1 &= (~i);
   69:             }
   70:         }
   71:         INSERT_WORDS(x,i0,i1);
   72:         return x;
   73: }
   74: 
   75: #if     LDBL_MANT_DIG == 53
   76: #ifdef  lint
   77: /* PROTOLIB1 */
   78: long double ceill(long double);
   79: #else   /* lint */
   80: __weak_alias(ceill, ceil);
   81: #endif  /* lint */
   82: #endif  /* LDBL_MANT_DIG == 53 */