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

diff --git a/src/port/rint.c b/src/port/rint.c
index 9c4d77535729384560bec83b0b231999205fffc1..d27fdfa6b4a1b89529f1962456cb0896f6d03575 100644 (file)
--- 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 <float.h>
 #include <math.h>
 
+/*
+ * 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 53b484f718c67b7fdcc082445724446a51ff7233..311fe730a5b4a01fbca9c5e3b08c248bbefb32ff 100644 (file)
--- 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 fcb14e3855e0298196028b334db48e501b87ad96..83fe022d7fe83cceccbe87b4e168ffb15e837f90 100644 (file)
--- 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 a1c70ed3e8d3a6c6bcc30587c59abdb3a1c893f4..a4f0cc23ec5fb73c4149026a426b0b8a2988d4e0 100644 (file)
--- 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 e79c3a8af95c8372ebd8526f26fc504097868b42..da8be51886c23c50b183cae24a574e368a07dae1 100644 (file)
--- 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 bacfbb24ac2bd6387c13deca72f7d0ada64eba4b..5e9774e922722cfc31235c1e0ad524ba4cede37c 100644 (file)
--- 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 1843a6d33bc86e774190f00d23cf90cb3248680e..d1881405252716f596e6f718cbb161a31db050b9 100644 (file)
--- 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 2f7f30c91d384f6ffdf21db9ef19a4d748979677..69723759fe3993f8562526fc90b551000d4d6420 100644 (file)
--- 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);