Fix for bug #15936: "round" differs on Windows to Unix

Both of our own implementations of rint(3) were inconsistent with the
most common behavior of rint() on those platforms that have it: round
to nearest, break ties by rounding to nearest even.

Fixed by leaving just one implementation of rint() in our source tree,
and changing its behavior to match the most common native
implementations on other platforms.

6 files changed, 113 insertions, 36 deletions

diff --git a/configure.in b/configure.in
index 9591b5bbc5a..8c2c4a9cf79 100644
--- a/configure.in
+++ b/configure.in
@@ -825,7 +825,7 @@ AC_TYPE_SIZE_T
 AC_HEADER_DIRENT
 AC_HEADER_STDC
 AC_HEADER_SYS_WAIT
-AC_CHECK_HEADERS(fcntl.h float.h floatingpoint.h ieeefp.h limits.h \
+AC_CHECK_HEADERS(fcntl.h fenv.h float.h floatingpoint.h ieeefp.h limits.h \
  memory.h pwd.h select.h \
  stdlib.h stddef.h \
  strings.h string.h synch.h sys/mman.h sys/socket.h netinet/in.h arpa/inet.h \
@@ -2060,7 +2060,7 @@ AC_FUNC_UTIME_NULL
 AC_FUNC_VPRINTF
 
 AC_CHECK_FUNCS(alarm bcmp bfill bmove bzero chsize cuserid fchmod fcntl \
-  fconvert fdatasync finite fpresetsticky fpsetmask fsync ftruncate \
+  fconvert fdatasync fesetround finite fpresetsticky fpsetmask fsync ftruncate \
   getcwd gethostbyaddr_r gethostbyname_r getpass getpassphrase getpwnam \
   getpwuid getrlimit getrusage getwd gmtime_r index initgroups isnan \
   localtime_r locking longjmp lrand48 madvise mallinfo memcpy memmove \
diff --git a/include/config-win.h b/include/config-win.h
index eba699159c7..ab463a7c142 100644
--- a/include/config-win.h
+++ b/include/config-win.h
@@ -31,7 +31,6 @@ functions */
 
 #include <sys/locking.h>
 #include <windows.h>
-#include <math.h>			/* Because of rint() */
 #include <fcntl.h>
 #include <io.h>
 #include <malloc.h>
@@ -223,13 +222,6 @@ typedef uint rf_SetTimer;
 #define inline __inline
 #endif /* __cplusplus */
 
-inline double rint(double nr)
-{
-  double f = floor(nr);
-  double c = ceil(nr);
-  return (((c-nr) >= (nr-f)) ? f :c);
-}
-
 #ifdef _WIN64
 #define ulonglong2double(A) ((double) (ulonglong) (A))
 #define my_off_t2double(A)  ((double) (my_off_t) (A))
@@ -281,7 +273,6 @@ inline ulonglong double2ulonglong(double d)
 #define HAVE_FLOAT_H
 #define HAVE_LIMITS_H
 #define HAVE_STDDEF_H
-#define HAVE_RINT		/* defined in this file */
 #define NO_FCNTL_NONBLOCK	/* No FCNTL */
 #define HAVE_ALLOCA
 #define HAVE_STRPBRK
diff --git a/include/my_global.h b/include/my_global.h
index 845ae042a42..3f872bfc855 100644
--- a/include/my_global.h
+++ b/include/my_global.h
@@ -484,8 +484,39 @@ typedef unsigned short ushort;
 #define set_bits(type, bit_count) (sizeof(type)*8 <= (bit_count) ? ~(type) 0 : ((((type) 1) << (bit_count)) - (type) 1))
 #define array_elements(A) ((uint) (sizeof(A)/sizeof(A[0])))
 #ifndef HAVE_RINT
-#define rint(A) floor((A)+(((A) < 0)? -0.5 : 0.5))
-#endif
+/**
+   All integers up to this number can be represented exactly as double precision
+   values (DBL_MANT_DIG == 53 for IEEE 754 hardware).
+*/
+#define MAX_EXACT_INTEGER ((1LL << DBL_MANT_DIG) - 1)
+
+/**
+   rint(3) implementation for platforms that do not have it.
+   Always rounds to the nearest integer with ties being rounded to the nearest
+   even integer to mimic glibc's rint() behavior in the "round-to-nearest"
+   FPU mode. Hardware-specific optimizations are possible (frndint on x86).
+   Unlike this implementation, hardware will also honor the FPU rounding mode.
+*/
+
+static inline double rint(double x)
+{
+  double f, i;
+  f = modf(x, &i);
+  /*
+    All doubles with absolute values > MAX_EXACT_INTEGER are even anyway,
+    no need to check it.
+  */
+  if (x > 0.0)
+    i += (double) ((f > 0.5) || (f == 0.5 &&
+                                 i <= (double) MAX_EXACT_INTEGER &&
+                                 (longlong) i % 2));
+  else
+    i -= (double) ((f < -0.5) || (f == -0.5 &&
+                                  i >= (double) -MAX_EXACT_INTEGER &&
+                                  (longlong) i % 2));
+  return i;
+}
+#endif /* HAVE_RINT */
 
 /* Define some general constants */
 #ifndef TRUE
diff --git a/mysql-test/r/func_math.result b/mysql-test/r/func_math.result
index 0d7adbbba5e..87cfb5b86a5 100644
--- a/mysql-test/r/func_math.result
+++ b/mysql-test/r/func_math.result
@@ -360,4 +360,34 @@ SELECT a DIV 2 FROM t1 UNION SELECT a DIV 2 FROM t1;
 a DIV 2
 0
 DROP TABLE t1;
+CREATE TABLE t1 (a DOUBLE);
+INSERT INTO t1 VALUES (-1.1), (1.1),
+(-1.5), (1.5),
+(-1.9), (1.9),
+(-2.1), (2.1),
+(-2.5), (2.5),
+(-2.9), (2.9),
+# Check numbers with absolute values > 2^53 - 1 
+# (see comments for MAX_EXACT_INTEGER)
+(-1e16 - 0.5), (1e16 + 0.5),
+(-1e16 - 1.5), (1e16 + 1.5);
+SELECT a, ROUND(a) FROM t1;
+a	ROUND(a)
+-1.1	-1
+1.1	1
+-1.5	-2
+1.5	2
+-1.9	-2
+1.9	2
+-2.1	-2
+2.1	2
+-2.5	-2
+2.5	2
+-2.9	-3
+2.9	3
+-1e+16	-10000000000000000
+1e+16	10000000000000000
+-1e+16	-10000000000000002
+1e+16	10000000000000002
+DROP TABLE t1;
 End of 5.0 tests
diff --git a/mysql-test/t/func_math.test b/mysql-test/t/func_math.test
index 9f12fdd696e..593cfe90c1b 100644
--- a/mysql-test/t/func_math.test
+++ b/mysql-test/t/func_math.test
@@ -229,5 +229,25 @@ INSERT INTO t1 VALUES ('a');
 SELECT a DIV 2 FROM t1 UNION SELECT a DIV 2 FROM t1;
 DROP TABLE t1;
 
+#
+# Bug #15936: "round" differs on Windows to Unix
+#
+
+CREATE TABLE t1 (a DOUBLE);
+
+INSERT INTO t1 VALUES (-1.1), (1.1),
+                      (-1.5), (1.5),
+                      (-1.9), (1.9),
+                      (-2.1), (2.1),
+                      (-2.5), (2.5),
+                      (-2.9), (2.9),
+# Check numbers with absolute values > 2^53 - 1 
+# (see comments for MAX_EXACT_INTEGER)
+                      (-1e16 - 0.5), (1e16 + 0.5),
+                      (-1e16 - 1.5), (1e16 + 1.5);
+
+SELECT a, ROUND(a) FROM t1;
+
+DROP TABLE t1;
 
 --echo End of 5.0 tests
diff --git a/sql/mysqld.cc b/sql/mysqld.cc
index ca68976d939..7856309b095 100644
--- a/sql/mysqld.cc
+++ b/sql/mysqld.cc
@@ -186,39 +186,44 @@ int initgroups(const char *,unsigned int);
 #ifdef HAVE_FP_EXCEPT				// Fix type conflict
 typedef fp_except fp_except_t;
 #endif
+#endif /* __FreeBSD__ && HAVE_IEEEFP_H */
+#ifdef HAVE_FENV_H
+#include <fenv.h>
+#endif
+#ifdef HAVE_SYS_FPU_H
+/* for IRIX to use set_fpc_csr() */
+#include <sys/fpu.h>
+#endif
 
+inline void setup_fpu()
+{
+#if defined(__FreeBSD__) && defined(HAVE_IEEEFP_H)
   /* We can't handle floating point exceptions with threads, so disable
      this on freebsd
+     Don't fall for overflow, underflow,divide-by-zero or loss of precision
   */
-
-inline void set_proper_floating_point_mode()
-{
-  /* Don't fall for overflow, underflow,divide-by-zero or loss of precision */
 #if defined(__i386__)
   fpsetmask(~(FP_X_INV | FP_X_DNML | FP_X_OFL | FP_X_UFL | FP_X_DZ |
 	      FP_X_IMP));
 #else
- fpsetmask(~(FP_X_INV |             FP_X_OFL | FP_X_UFL | FP_X_DZ |
-	     FP_X_IMP));
-#endif
-}
-#elif defined(__sgi)
-/* for IRIX to use set_fpc_csr() */
-#include <sys/fpu.h>
+  fpsetmask(~(FP_X_INV |             FP_X_OFL | FP_X_UFL | FP_X_DZ |
+              FP_X_IMP));
+#endif /* __i386__ */
+#endif /* __FreeBSD__ && HAVE_IEEEFP_H */
 
-inline void set_proper_floating_point_mode()
-{
+#ifdef HAVE_FESETROUND
+    /* Set FPU rounding mode to "round-to-nearest" */
+  fesetround(FE_TONEAREST);
+#endif /* HAVE_FESETROUND */
+    
+#if defined(__sgi) && defined(HAVE_SYS_FPU_H)
   /* Enable denormalized DOUBLE values support for IRIX */
-  {
-    union fpc_csr n;
-    n.fc_word = get_fpc_csr();
-    n.fc_struct.flush = 0;
-    set_fpc_csr(n.fc_word);
-  }
+  union fpc_csr n;
+  n.fc_word = get_fpc_csr();
+  n.fc_struct.flush = 0;
+  set_fpc_csr(n.fc_word);
+#endif
 }
-#else
-#define set_proper_floating_point_mode()
-#endif /* __FreeBSD__ && HAVE_IEEEFP_H */
 
 } /* cplusplus */
 
@@ -3279,7 +3284,7 @@ static int init_server_components()
   query_cache_init();
   query_cache_resize(query_cache_size);
   randominit(&sql_rand,(ulong) server_start_time,(ulong) server_start_time/2);
-  set_proper_floating_point_mode();
+  setup_fpu();
   init_thr_lock();
 #ifdef HAVE_REPLICATION
   init_slave_list();