[BACK]Return to s_rint.c CVS log [TXT][DIR] Up to [cvs.NetBSD.org] / src / lib / libm / src

Annotation of src/lib/libm/src/s_rint.c, Revision 1.6

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

CVSweb <webmaster@jp.NetBSD.org>