gonzui

    1: /* @(#)s_rint.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: /*
   14:  * rint(x)
   15:  * Return x rounded to integral value according to the prevailing
   16:  * rounding mode.
   17:  * Method:
   18:  *      Using floating addition.
   19:  * Exception:
   20:  *      Inexact flag raised if x not equal to rint(x).
   21:  */
   22: 
   23: /* LINTLIBRARY */
   24: 
   25: #include <sys/cdefs.h>
   26: #include <float.h>
   27: #include <math.h>
   28: 
   29: #include "math_private.h"
   30: 
   31: static const double
   32: TWO52[2]={
   33:   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
   34:  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
   35: };
   36: 
   37: double
   38: rint(double x)
   39: {
   40:         int32_t i0,jj0,sx;
   41:         u_int32_t i,i1;
   42:         double t;
   43:         volatile double w;     /* clip extra precision */
   44:         EXTRACT_WORDS(i0,i1,x);
   45:         sx = (i0>>31)&1;
   46:         jj0 = ((i0>>20)&0x7ff)-0x3ff;
   47:         if(jj0<20) {
   48:             if(jj0<0) {        
   49:                 if(((i0&0x7fffffff)|i1)==0) return x;
   50:                 i1 |= (i0&0x0fffff);
   51:                 i0 &= 0xfffe0000;
   52:                 i0 |= ((i1|-i1)>>12)&0x80000;
   53:                 SET_HIGH_WORD(x,i0);
   54:                 w = TWO52[sx]+x;
   55:                 t =  w-TWO52[sx];
   56:                 GET_HIGH_WORD(i0,t);
   57:                 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
   58:                 return t;
   59:             } else {
   60:                 i = (0x000fffff)>>jj0;
   61:                 if(((i0&i)|i1)==0) return x; /* x is integral */
   62:                 i>>=1;
   63:                 if(((i0&i)|i1)!=0) {
   64:                     if(jj0==19) i1 = 0x40000000; else
   65:                     i0 = (i0&(~i))|((0x20000)>>jj0);
   66:                 }
   67:             }
   68:         } else if (jj0>51) {
   69:             if(jj0==0x400) return x+x; /* inf or NaN */
   70:             else return x;             /* x is integral */
   71:         } else {
   72:             i = ((u_int32_t)(0xffffffff))>>(jj0-20);
   73:             if((i1&i)==0) return x;    /* x is integral */
   74:             i>>=1;
   75:             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(jj0-20));
   76:         }
   77:         INSERT_WORDS(x,i0,i1);
   78:         w = TWO52[sx]+x;
   79:         return w-TWO52[sx];
   80: }
   81: 
   82: #if     LDBL_MANT_DIG == 53
   83: #ifdef  lint
   84: /* PROTOLIB1 */
   85: long double rintl(long double);
   86: #else   /* lint */
   87: __weak_alias(rintl, rint);
   88: #endif  /* lint */
   89: #endif  /* LDBL_MANT_DIG == 53 */