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