|
3 | 3 | * rint.c |
4 | 4 | * rint() implementation |
5 | 5 | * |
| 6 | + * By Pedro Gimeno Fortea, donated to the public domain |
| 7 | + * |
6 | 8 | * IDENTIFICATION |
7 | 9 | * src/port/rint.c |
8 | 10 | * |
9 | 11 | *------------------------------------------------------------------------- |
10 | 12 | */ |
11 | 13 | #include"c.h" |
12 | 14 |
|
| 15 | +#include<float.h> |
13 | 16 | #include<math.h> |
14 | 17 |
|
| 18 | +/* |
| 19 | + * Round to nearest integer, with halfway cases going to the nearest even. |
| 20 | + */ |
15 | 21 | double |
16 | 22 | rint(doublex) |
17 | 23 | { |
18 | | -return (x >=0.0) ?floor(x+0.5) :ceil(x-0.5); |
| 24 | +doublex_orig; |
| 25 | +doubler; |
| 26 | + |
| 27 | +/* Per POSIX, NaNs must be returned unchanged. */ |
| 28 | +if (isnan(x)) |
| 29 | +returnx; |
| 30 | + |
| 31 | +if (x <=0.0) |
| 32 | +{ |
| 33 | +/* Both positive and negative zero should be returned unchanged. */ |
| 34 | +if (x==0.0) |
| 35 | +returnx; |
| 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 | +returnx_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 | +returnr; |
| 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 | +returnfloor(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 | +return0.0; |
| 91 | +if (x==x_orig-1.0) |
| 92 | +returnx_orig; |
| 93 | +r=ceil(x); |
| 94 | +if (r!=x) |
| 95 | +returnr; |
| 96 | +returnceil(x*0.5)*2.0; |
| 97 | +} |
19 | 98 | } |