1 /*
2  * Copyright (c) 1998, 2001, Oracle and/or its affiliates. All rights reserved.
3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4  *
5  * This code is free software; you can redistribute it and/or modify it
6  * under the terms of the GNU General Public License version 2 only, as
7  * published by the Free Software Foundation.  Oracle designates this
8  * particular file as subject to the "Classpath" exception as provided
9  * by Oracle in the LICENSE file that accompanied this code.
10  *
11  * This code is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14  * version 2 for more details (a copy is included in the LICENSE file that
15  * accompanied this code).
16  *
17  * You should have received a copy of the GNU General Public License version
18  * 2 along with this work; if not, write to the Free Software Foundation,
19  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20  *
21  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22  * or visit www.oracle.com if you need additional information or have any
23  * questions.
24  */
25 
26 /*
27  * rint(x)
28  * Return x rounded to integral value according to the prevailing
29  * rounding mode.
30  * Method:
31  *      Using floating addition.
32  * Exception:
33  *      Inexact flag raised if x not equal to rint(x).
34  */
35 
36 #include "fdlibm.h"
37 
38 #ifdef __STDC__
39 static const double
40 #else
41 static double
42 #endif
43 TWO52[2]={
44   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
45  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
46 };
47 
48 #ifdef __STDC__
rint(double x)49         double rint(double x)
50 #else
51         double rint(x)
52         double x;
53 #endif
54 {
55         int i0,j0,sx;
56         unsigned i,i1;
57         double w,t;
58         i0 =  __HI(x);
59         sx = (i0>>31)&1;
60         i1 =  __LO(x);
61         j0 = ((i0>>20)&0x7ff)-0x3ff;
62         if(j0<20) {
63             if(j0<0) {
64                 if(((i0&0x7fffffff)|i1)==0) return x;
65                 i1 |= (i0&0x0fffff);
66                 i0 &= 0xfffe0000;
67                 i0 |= ((i1|-i1)>>12)&0x80000;
68                 __HI(x)=i0;
69                 w = TWO52[sx]+x;
70                 t =  w-TWO52[sx];
71                 i0 = __HI(t);
72                 __HI(t) = (i0&0x7fffffff)|(sx<<31);
73                 return t;
74             } else {
75                 i = (0x000fffff)>>j0;
76                 if(((i0&i)|i1)==0) return x; /* x is integral */
77                 i>>=1;
78                 if(((i0&i)|i1)!=0) {
79                     if(j0==19) i1 = 0x40000000; else
80                     i0 = (i0&(~i))|((0x20000)>>j0);
81                 }
82             }
83         } else if (j0>51) {
84             if(j0==0x400) return x+x;   /* inf or NaN */
85             else return x;              /* x is integral */
86         } else {
87             i = ((unsigned)(0xffffffff))>>(j0-20);
88             if((i1&i)==0) return x;     /* x is integral */
89             i>>=1;
90             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
91         }
92         __HI(x) = i0;
93         __LO(x) = i1;
94         w = TWO52[sx]+x;
95         return w-TWO52[sx];
96 }
97