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>