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>