1 /*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12 /*
13 * rint(x)
14 * Return x rounded to integral value according to the prevailing
15 * rounding mode.
16 * Method:
17 * Using floating addition.
18 * Exception:
19 * Inexact flag raised if x not equal to rint(x).
20 */
21
22 #include <float.h>
23
24 #include "math.h"
25 #include "math_private.h"
26
27 static const double
28 TWO52[2]={
29 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
30 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
31 };
32
33 double
rint(double x)34 rint(double x)
35 {
36 int32_t i0,j0,sx;
37 u_int32_t i,i1;
38 double w,t;
39 EXTRACT_WORDS(i0,i1,x);
40 sx = (i0>>31)&1;
41 j0 = ((i0>>20)&0x7ff)-0x3ff;
42 if(j0<20) {
43 if(j0<0) {
44 if(((i0&0x7fffffff)|i1)==0) return x;
45 i1 |= (i0&0x0fffff);
46 i0 &= 0xfffe0000;
47 i0 |= ((i1|-i1)>>12)&0x80000;
48 SET_HIGH_WORD(x,i0);
49 STRICT_ASSIGN(double,w,TWO52[sx]+x);
50 t = w-TWO52[sx];
51 GET_HIGH_WORD(i0,t);
52 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
53 return t;
54 } else {
55 i = (0x000fffff)>>j0;
56 if(((i0&i)|i1)==0) return x; /* x is integral */
57 i>>=1;
58 if(((i0&i)|i1)!=0) {
59 /*
60 * Some bit is set after the 0.5 bit. To avoid the
61 * possibility of errors from double rounding in
62 * w = TWO52[sx]+x, adjust the 0.25 bit to a lower
63 * guard bit. We do this for all j0<=51. The
64 * adjustment is trickiest for j0==18 and j0==19
65 * since then it spans the word boundary.
66 */
67 if(j0==19) i1 = 0x40000000; else
68 if(j0==18) i1 = 0x80000000; else
69 i0 = (i0&(~i))|((0x20000)>>j0);
70 }
71 }
72 } else if (j0>51) {
73 if(j0==0x400) return x+x; /* inf or NaN */
74 else return x; /* x is integral */
75 } else {
76 i = ((u_int32_t)(0xffffffff))>>(j0-20);
77 if((i1&i)==0) return x; /* x is integral */
78 i>>=1;
79 if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
80 }
81 INSERT_WORDS(x,i0,i1);
82 STRICT_ASSIGN(double,w,TWO52[sx]+x);
83 return w-TWO52[sx];
84 }
85
86 #if (LDBL_MANT_DIG == 53)
87 __weak_reference(rint, rintl);
88 #endif
89