1 /*-------------------------------------------------------------------------
2  *
3  * rint.c
4  *	  rint() implementation
5  *
6  * By Pedro Gimeno Fortea, donated to the public domain
7  *
8  * IDENTIFICATION
9  *	  src/port/rint.c
10  *
11  *-------------------------------------------------------------------------
12  */
13 #include "c.h"
14 
15 #include <float.h>
16 #include <math.h>
17 
18 /*
19  * Round to nearest integer, with halfway cases going to the nearest even.
20  */
21 double
22 rint(double x)
23 {
24 	double		x_orig;
25 	double		r;
26 
27 	/* Per POSIX, NaNs must be returned unchanged. */
28 	if (isnan(x))
29 		return x;
30 
31 	if (x <= 0.0)
32 	{
33 		/* Both positive and negative zero should be returned unchanged. */
34 		if (x == 0.0)
35 			return x;
36 
37 		/*
38 		 * Subtracting 0.5 from a number very close to -0.5 can round to
39 		 * exactly -1.0, producing incorrect results, so we take the opposite
40 		 * approach: add 0.5 to the negative number, so that it goes closer to
41 		 * zero (or at most to +0.5, which is dealt with next), avoiding the
42 		 * precision issue.
43 		 */
44 		x_orig = x;
45 		x += 0.5;
46 
47 		/*
48 		 * Be careful to return minus zero when input+0.5 >= 0, as that's what
49 		 * rint() should return with negative input.
50 		 */
51 		if (x >= 0.0)
52 			return -0.0;
53 
54 		/*
55 		 * For very big numbers the input may have no decimals.  That case is
56 		 * detected by testing x+0.5 == x+1.0; if that happens, the input is
57 		 * returned unchanged.  This also covers the case of minus infinity.
58 		 */
59 		if (x == x_orig + 1.0)
60 			return x_orig;
61 
62 		/* Otherwise produce a rounded estimate. */
63 		r = floor(x);
64 
65 		/*
66 		 * If the rounding did not produce exactly input+0.5 then we're done.
67 		 */
68 		if (r != x)
69 			return r;
70 
71 		/*
72 		 * The original fractional part was exactly 0.5 (since
73 		 * floor(input+0.5) == input+0.5).  We need to round to nearest even.
74 		 * Dividing input+0.5 by 2, taking the floor and multiplying by 2
75 		 * yields the closest even number.  This part assumes that division by
76 		 * 2 is exact, which should be OK because underflow is impossible
77 		 * here: x is an integer.
78 		 */
79 		return floor(x * 0.5) * 2.0;
80 	}
81 	else
82 	{
83 		/*
84 		 * The positive case is similar but with signs inverted and using
85 		 * ceil() instead of floor().
86 		 */
87 		x_orig = x;
88 		x -= 0.5;
89 		if (x <= 0.0)
90 			return 0.0;
91 		if (x == x_orig - 1.0)
92 			return x_orig;
93 		r = ceil(x);
94 		if (r != x)
95 			return r;
96 		return ceil(x * 0.5) * 2.0;
97 	}
98 }
99