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
27static const double
28TWO52[2]={
29  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
30 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
31};
32
33double
34rint(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