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

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.10    ! simonb     15: __RCSID("$NetBSD: s_rint.c,v 1.9 1997/10/09 11:33:10 lukem 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: #ifdef __STDC__
                     32: static const double
                     33: #else
1.10    ! simonb     34: static double
1.1       jtc        35: #endif
                     36: TWO52[2]={
                     37:   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
                     38:  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
                     39: };
                     40:
                     41: #ifdef __STDC__
                     42:        double rint(double x)
                     43: #else
                     44:        double rint(x)
                     45:        double x;
                     46: #endif
                     47: {
1.6       jtc        48:        int32_t i0,j0,sx;
                     49:        u_int32_t i,i1;
1.1       jtc        50:        double w,t;
1.5       jtc        51:        EXTRACT_WORDS(i0,i1,x);
1.1       jtc        52:        sx = (i0>>31)&1;
                     53:        j0 = ((i0>>20)&0x7ff)-0x3ff;
                     54:        if(j0<20) {
1.10    ! simonb     55:            if(j0<0) {
1.1       jtc        56:                if(((i0&0x7fffffff)|i1)==0) return x;
                     57:                i1 |= (i0&0x0fffff);
                     58:                i0 &= 0xfffe0000;
                     59:                i0 |= ((i1|-i1)>>12)&0x80000;
1.5       jtc        60:                SET_HIGH_WORD(x,i0);
1.1       jtc        61:                w = TWO52[sx]+x;
                     62:                t =  w-TWO52[sx];
1.5       jtc        63:                GET_HIGH_WORD(i0,t);
                     64:                SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
1.1       jtc        65:                return t;
                     66:            } else {
                     67:                i = (0x000fffff)>>j0;
                     68:                if(((i0&i)|i1)==0) return x; /* x is integral */
                     69:                i>>=1;
                     70:                if(((i0&i)|i1)!=0) {
                     71:                    if(j0==19) i1 = 0x40000000; else
                     72:                    i0 = (i0&(~i))|((0x20000)>>j0);
                     73:                }
                     74:            }
                     75:        } else if (j0>51) {
                     76:            if(j0==0x400) return x+x;   /* inf or NaN */
                     77:            else return x;              /* x is integral */
                     78:        } else {
1.6       jtc        79:            i = ((u_int32_t)(0xffffffff))>>(j0-20);
1.1       jtc        80:            if((i1&i)==0) return x;     /* x is integral */
                     81:            i>>=1;
                     82:            if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
                     83:        }
1.5       jtc        84:        INSERT_WORDS(x,i0,i1);
1.1       jtc        85:        w = TWO52[sx]+x;
                     86:        return w-TWO52[sx];
                     87: }

CVSweb <webmaster@jp.NetBSD.org>