From 06bf0dd6e354c765403d1331cc9896b360754521 Mon Sep 17 00:00:00 2001 From: Tom Lane Date: Wed, 25 Mar 2015 15:54:08 -0400 Subject: [PATCH] Upgrade src/port/rint.c to be POSIX-compliant. The POSIX spec says that rint() rounds halfway cases to nearest even. Our substitute implementation failed to do that, rather rounding halfway cases away from zero; and it also got some other cases (such as minus zero) wrong. This led to observable cross-platform differences, as reported in bug #12885 from Rich Schaaf; in particular, casting from float to int didn't honor round-to-nearest-even on builds using rint.c. Implement something that attempts to cover all cases per spec, and add some simple regression tests so that we'll notice if any platforms still get this wrong. Although this is a bug fix, no back-patch, as a behavioral change in the back branches was agreed not to be a good idea. Pedro Gimeno Fortea, reviewed by Michael Paquier and myself --- src/port/rint.c | 81 ++++++++++++++++++- src/test/regress/expected/int2.out | 20 +++++ src/test/regress/expected/int4.out | 20 +++++ .../expected/int8-exp-three-digits.out | 20 +++++ src/test/regress/expected/int8.out | 20 +++++ src/test/regress/sql/int2.sql | 10 +++ src/test/regress/sql/int4.sql | 10 +++ src/test/regress/sql/int8.sql | 10 +++ 8 files changed, 190 insertions(+), 1 deletion(-) diff --git a/src/port/rint.c b/src/port/rint.c index 9c4d775357..d27fdfa6b4 100644 --- a/src/port/rint.c +++ b/src/port/rint.c @@ -3,6 +3,8 @@ * rint.c * rint() implementation * + * By Pedro Gimeno Fortea, donated to the public domain + * * IDENTIFICATION * src/port/rint.c * @@ -10,10 +12,87 @@ */ #include "c.h" +#include #include +/* + * Round to nearest integer, with halfway cases going to the nearest even. + */ double rint(double x) { - return (x >= 0.0) ? floor(x + 0.5) : ceil(x - 0.5); + 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; + } } diff --git a/src/test/regress/expected/int2.out b/src/test/regress/expected/int2.out index 53b484f718..311fe730a5 100644 --- a/src/test/regress/expected/int2.out +++ b/src/test/regress/expected/int2.out @@ -266,3 +266,23 @@ SELECT (-32768)::int2 % (-1)::int2; 0 (1 row) +-- check rounding when casting from float +SELECT x, x::int2 AS int2_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x); + x | int2_value +------+------------ + -2.5 | -2 + -1.5 | -2 + -0.5 | 0 + 0 | 0 + 0.5 | 0 + 1.5 | 2 + 2.5 | 2 +(7 rows) + diff --git a/src/test/regress/expected/int4.out b/src/test/regress/expected/int4.out index fcb14e3855..83fe022d7f 100644 --- a/src/test/regress/expected/int4.out +++ b/src/test/regress/expected/int4.out @@ -363,3 +363,23 @@ SELECT (-2147483648)::int4 % (-1)::int2; 0 (1 row) +-- check rounding when casting from float +SELECT x, x::int4 AS int4_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x); + x | int4_value +------+------------ + -2.5 | -2 + -1.5 | -2 + -0.5 | 0 + 0 | 0 + 0.5 | 0 + 1.5 | 2 + 2.5 | 2 +(7 rows) + diff --git a/src/test/regress/expected/int8-exp-three-digits.out b/src/test/regress/expected/int8-exp-three-digits.out index a1c70ed3e8..a4f0cc23ec 100644 --- a/src/test/regress/expected/int8-exp-three-digits.out +++ b/src/test/regress/expected/int8-exp-three-digits.out @@ -846,3 +846,23 @@ SELECT (-9223372036854775808)::int8 % (-1)::int2; 0 (1 row) +-- check rounding when casting from float +SELECT x, x::int8 AS int8_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x); + x | int8_value +------+------------ + -2.5 | -2 + -1.5 | -2 + -0.5 | 0 + 0 | 0 + 0.5 | 0 + 1.5 | 2 + 2.5 | 2 +(7 rows) + diff --git a/src/test/regress/expected/int8.out b/src/test/regress/expected/int8.out index e79c3a8af9..da8be51886 100644 --- a/src/test/regress/expected/int8.out +++ b/src/test/regress/expected/int8.out @@ -846,3 +846,23 @@ SELECT (-9223372036854775808)::int8 % (-1)::int2; 0 (1 row) +-- check rounding when casting from float +SELECT x, x::int8 AS int8_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x); + x | int8_value +------+------------ + -2.5 | -2 + -1.5 | -2 + -0.5 | 0 + 0 | 0 + 0.5 | 0 + 1.5 | 2 + 2.5 | 2 +(7 rows) + diff --git a/src/test/regress/sql/int2.sql b/src/test/regress/sql/int2.sql index bacfbb24ac..5e9774e922 100644 --- a/src/test/regress/sql/int2.sql +++ b/src/test/regress/sql/int2.sql @@ -92,3 +92,13 @@ SELECT ((-1::int2<<15)+1::int2)::text; SELECT (-32768)::int2 * (-1)::int2; SELECT (-32768)::int2 / (-1)::int2; SELECT (-32768)::int2 % (-1)::int2; + +-- check rounding when casting from float +SELECT x, x::int2 AS int2_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x); diff --git a/src/test/regress/sql/int4.sql b/src/test/regress/sql/int4.sql index 1843a6d33b..d188140525 100644 --- a/src/test/regress/sql/int4.sql +++ b/src/test/regress/sql/int4.sql @@ -135,3 +135,13 @@ SELECT (-2147483648)::int4 % (-1)::int4; SELECT (-2147483648)::int4 * (-1)::int2; SELECT (-2147483648)::int4 / (-1)::int2; SELECT (-2147483648)::int4 % (-1)::int2; + +-- check rounding when casting from float +SELECT x, x::int4 AS int4_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x); diff --git a/src/test/regress/sql/int8.sql b/src/test/regress/sql/int8.sql index 2f7f30c91d..69723759fe 100644 --- a/src/test/regress/sql/int8.sql +++ b/src/test/regress/sql/int8.sql @@ -205,3 +205,13 @@ SELECT (-9223372036854775808)::int8 % (-1)::int4; SELECT (-9223372036854775808)::int8 * (-1)::int2; SELECT (-9223372036854775808)::int8 / (-1)::int2; SELECT (-9223372036854775808)::int8 % (-1)::int2; + +-- check rounding when casting from float +SELECT x, x::int8 AS int8_value +FROM (VALUES (-2.5::float8), + (-1.5::float8), + (-0.5::float8), + (0.0::float8), + (0.5::float8), + (1.5::float8), + (2.5::float8)) t(x);