[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.11

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

CVSweb <webmaster@jp.NetBSD.org>