gonzui

    1: /* @(#)s_floor.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:  * floor(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 floor(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: floor(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=i1=0;} 
   43:                     else if(((i0&0x7fffffff)|i1)!=0)
   44:                         { i0=0xbff00000;i1=0;}
   45:                 }
   46:             } else {
   47:                 i = (0x000fffff)>>jj0;
   48:                 if(((i0&i)|i1)==0) return x; /* x is integral */
   49:                 if(huge+x>0.0) {      /* raise inexact flag */
   50:                     if(i0<0) i0 += (0x00100000)>>jj0;
   51:                     i0 &= (~i); i1=0;
   52:                 }
   53:             }
   54:         } else if (jj0>51) {
   55:             if(jj0==0x400) return x+x; /* inf or NaN */
   56:             else return x;             /* x is integral */
   57:         } else {
   58:             i = ((u_int32_t)(0xffffffff))>>(jj0-20);
   59:             if((i1&i)==0) return x;    /* x is integral */
   60:             if(huge+x>0.0) {           /* raise inexact flag */
   61:                 if(i0<0) {
   62:                     if(jj0==20) i0+=1; 
   63:                     else {
   64:                         j = i1+(1<<(52-jj0));
   65:                         if(j<i1) i0 +=1 ;    /* got a carry */
   66:                         i1=j;
   67:                     }
   68:                 }
   69:                 i1 &= (~i);
   70:             }
   71:         }
   72:         INSERT_WORDS(x,i0,i1);
   73:         return x;
   74: }
   75: 
   76: #if     LDBL_MANT_DIG == 53
   77: #ifdef  lint
   78: /* PROTOLIB1 */
   79: long double floorl(long double);
   80: #else   /* lint */
   81: __weak_alias(floorl, floor);
   82: #endif  /* lint */
   83: #endif  /* LDBL_MANT_DIG == 53 */