Avoid precision errors in round(), lround() and llround() functions.

14 files changed, 325 insertions, 150 deletions

diff --git a/winsup/mingw/ChangeLog b/winsup/mingw/ChangeLog
index e64f943bc..0cc75f7df 100644
--- a/winsup/mingw/ChangeLog
+++ b/winsup/mingw/ChangeLog
@@ -1,3 +1,34 @@
+2008-06-02  Keith Marshall  <keithmarshall@users.sourceforge.net>
+
+	Avoid precision errors in round(), lround() and llround() functions.
+
+	* mingwex/math/round_generic.c: New file; it replaces...
+	* mingwex/math/round.c: ...this; delete it.
+	* mingwex/math/roundf.c: ...and this; ditto.
+	* mingwex/math/roundl.c: ...and this; ditto.
+
+	* mingwex/math/lround_generic.c: New file; it replaces...
+	* mingwex/math/lround.c: ...this; delete it.
+	* mingwex/math/lroundf.c: ...and this; ditto.
+	* mingwex/math/lroundl.c: ...and this; ditto.
+	* mingwex/math/llround.c: ...and this; ditto.
+	* mingwex/math/llroundf.c: ...and this; ditto.
+	* mingwex/math/llroundl.c: ...and this; ditto.
+
+	* mingwex/math/round_internal.h: New file; it provides generic common
+	code, which is shared by `round_generic.c' and `lround_generic.c'; the
+	implementation is based on techniques suggested by Danny Smith and
+	Greg Chicares.
+
+	* mingwex/Makefile.in (MATH_DISTFILES): Remove `round.c', `roundf.c',
+	`roundl.c', `lround.c', `lroundf.c', `lroundl.c', `llround.c',
+	`llroundf.c' and `llroundl.c'; replace by `round_internal.h',
+	`round_generic.c' and `lround_generic.c'.
+	(MATH_OBJS): Factor out files listed in...
+	(MATH_ROUND_OBJS, MATH_LROUND_OBJS, MATH_LLROUND_OBJS): ...these new
+	macros; define them; specify dependencies and build rules; add to...
+	(LIB_OBJS): ...this list.
+
 2008-05-22  Danny Smith  <dannysmith@users.sourceforge.net>
 
 	[ mingw-Bugs-1961893 ]
diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in
index 9010634cb..9563ad4e7 100644
--- a/winsup/mingw/mingwex/Makefile.in
+++ b/winsup/mingw/mingwex/Makefile.in
@@ -54,17 +54,16 @@ MATH_DISTFILES = \
 	fpclassify.c fpclassifyf.c fpclassifyl.c \
 	frexpf.c frexpl.S fucom.c hypotf.c hypotl.c ilogb.S ilogbf.S \
 	ilogbl.S isnan.c isnanf.c isnanl.c ldexpf.c ldexpl.c \
-	lgamma.c lgammaf.c lgammal.c llrint.c \
-	llrintf.c llrintl.c llround.c llroundf.c llroundl.c \
+	lgamma.c lgammaf.c lgammal.c llrint.c llrintf.c llrintl.c \
 	log10f.S log10l.S log1p.S log1pf.S log1pl.S log2.S log2f.S \
 	log2l.S logb.c logbf.c logbl.c logf.S logl.S lrint.c lrintf.c \
-	lrintl.c lround.c lroundf.c lroundl.c modff.c modfl.c \
+	lrintl.c lround_generic.c modff.c modfl.c \
 	nearbyint.S nearbyintf.S nearbyintl.S \
 	nextafterf.c nextafterl.c nexttowardf.c nexttoward.c \
 	powf.c powi.c powif.c powil.c powl.c \
 	remainder.S remainderf.S remainderl.S remquo.S \
-	remquof.S remquol.S rint.c rintf.c rintl.c round.c roundf.c \
-	roundl.c scalbn.S scalbnf.S scalbnl.S s_erf.c sf_erf.c \
+	remquof.S remquol.S rint.c rintf.c rintl.c round_internal.h \
+	round_generic.c scalbn.S scalbnf.S scalbnl.S s_erf.c sf_erf.c \
 	signbit.c signbitf.c signbitl.c sinf.S sinhf.c sinhl.c sinl.S \
 	sqrtf.c sqrtl.c tanf.S tanhf.c tanhl.c tanl.S tgamma.c \
 	tgammaf.c tgammal.c trunc.c truncf.c truncl.c \
@@ -153,22 +152,24 @@ MATH_OBJS = \
 	fpclassify.o fpclassifyf.o fpclassifyl.o \
 	frexpf.o frexpl.o fucom.o hypotf.o hypotl.o ilogb.o ilogbf.o \
 	ilogbl.o isnan.o isnanf.o isnanl.o ldexpf.o ldexpl.o \
-	lgamma.o lgammaf.o lgammal.o llrint.o \
-	llrintf.o llrintl.o llround.o llroundf.o llroundl.o \
+	lgamma.o lgammaf.o lgammal.o llrint.o llrintf.o llrintl.o \
 	log10f.o log10l.o log1p.o log1pf.o log1pl.o log2.o log2f.o \
-	log2l.o logb.o logbf.o logbl.o logf.o logl.o lrint.o lrintf.o \
-	lrintl.o lround.o lroundf.o lroundl.o modff.o modfl.o \
+	log2l.o logb.o logbf.o logbl.o logf.o logl.o \
+	lrint.o lrintf.o lrintl.o modff.o modfl.o \
 	nearbyint.o nearbyintf.o nearbyintl.o \
 	nextafterf.o nextafterl.o nexttowardf.o nexttoward.o \
 	powf.o powi.o powif.o powil.o powl.o \
 	remainder.o remainderf.o remainderl.o remquo.o \
-	remquof.o remquol.o rint.o rintf.o rintl.o round.o roundf.o \
-	roundl.o scalbn.o scalbnf.o scalbnl.o s_erf.o sf_erf.o \
+	remquof.o remquol.o rint.o rintf.o rintl.o \
+	scalbn.o scalbnf.o scalbnl.o s_erf.o sf_erf.o \
 	signbit.o signbitf.o signbitl.o sinf.o sinhf.o sinhl.o sinl.o \
 	sqrtf.o sqrtl.o tanf.o tanhf.o tanhl.o tanl.o tgamma.o \
 	tgammaf.o tgammal.o trunc.o truncf.o truncl.o \
 	acosh.o acoshf.o acoshl.o asinh.o asinhf.o asinhl.o \
 	atanh.o atanhf.o atanhl.o
+MATH_ROUND_OBJS = round.o roundf.o roundl.o
+MATH_LROUND_OBJS = lround.o lroundf.o lroundl.o
+MATH_LLROUND_OBJS = llround.o llroundf.o llroundl.o
 FENV_OBJS = fesetround.o  fegetround.o \
 	fegetenv.o fesetenv.o feupdateenv.o \
 	feclearexcept.o feholdexcept.o fegetexceptflag.o \
@@ -196,7 +197,8 @@ GDTOA_OBJS = \
 	mingw_snprintf.o
 
 LIB_OBJS = $(Q8_OBJS)  $(CTYPE_OBJS) $(STDLIB_STUB_OBJS) \
-	$(STDIO_OBJS) $(MATH_OBJS)  $(FENV_OBJS) \
+	$(STDIO_OBJS) $(MATH_OBJS) $(MATH_ROUND_OBJS) \
+	$(MATH_LROUND_OBJS) $(MATH_LLROUND_OBJS) $(FENV_OBJS) \
 	$(POSIX_OBJS) $(REPLACE_OBJS) $(COMPLEX_OBJS) \
 	$(GDTOA_OBJS)
 
@@ -210,6 +212,14 @@ $(LIBMINGWEX_A): $(LIB_OBJS)
 	$(AR) $(ARFLAGS) $@ $(LIB_OBJS)
 	$(RANLIB) $@
 
+$(MATH_ROUND_OBJS): round_generic.c
+	$(CC) $(ALL_CFLAGS) -I$(srcdir)/math -c -o $@ \
+	  -D FUNCTION=$* $(srcdir)/math/round_generic.c
+
+$(MATH_LROUND_OBJS) $(MATH_LLROUND_OBJS): lround_generic.c
+	$(CC) $(ALL_CFLAGS) -I$(srcdir)/math -c -o $@ \
+	  -D FUNCTION=$* $(srcdir)/math/lround_generic.c
+
 Makefile: Makefile.in config.status configure
 	$(SHELL) config.status
 
@@ -255,6 +265,8 @@ mbrtowc.o wcrtomb.o wcstof.o wcstold.o: mb_wc_common.h
 
 $(GDTOA_OBJS): gd_arith.h gdtoa.h gdtoaimp.h gd_qnan.h
 
+$(MATH_ROUND_OBJS) $(MATH_LROUND_OBJS) $(MATH_LLROUND_OBJS): round_internal.h
+
 dist:
 	mkdir $(distdir)/mingwex
 	chmod 755 $(distdir)/mingwex
diff --git a/winsup/mingw/mingwex/math/llround.c b/winsup/mingw/mingwex/math/llround.c
deleted file mode 100644
index 45b754c75..000000000
--- a/winsup/mingw/mingwex/math/llround.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <math.h>
-#include <limits.h>
-#include <errno.h>
-
-long long
-llround (double x)
-{
-  /* Add +/- 0.5, then round towards zero.  */
-  double tmp = trunc (x + (x >= 0.0 ?  0.5 : -0.5));
-  if (!isfinite (tmp) 
-      || tmp > (double)LONG_LONG_MAX
-      || tmp < (double)LONG_LONG_MIN)
-    { 
-      errno = ERANGE;
-      /* Undefined behaviour, so we could return anything.  */
-      /* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN;  */
-    }
-  return (long long)tmp;
-}
diff --git a/winsup/mingw/mingwex/math/llroundf.c b/winsup/mingw/mingwex/math/llroundf.c
deleted file mode 100644
index 6a6e9b51e..000000000
--- a/winsup/mingw/mingwex/math/llroundf.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <math.h>
-#include <limits.h>
-#include <errno.h>
-
-long long
-llroundf (float x)
-{
-  /* Add +/- 0.5, then round towards zero.  */
-  float tmp = truncf (x + (x >= 0.0F ?  0.5F : -0.5F));
-  if (!isfinite (tmp) 
-      || tmp > (float)LONG_LONG_MAX
-      || tmp < (float)LONG_LONG_MIN)
-    { 
-      errno = ERANGE;
-      /* Undefined behaviour, so we could return anything.  */
-      /* return tmp > 0.0F ? LONG_LONG_MAX : LONG_LONG_MIN; */
-    }
-  return (long long)tmp;
-}  
diff --git a/winsup/mingw/mingwex/math/llroundl.c b/winsup/mingw/mingwex/math/llroundl.c
deleted file mode 100644
index 9d2217411..000000000
--- a/winsup/mingw/mingwex/math/llroundl.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <math.h>
-#include <limits.h>
-#include <errno.h>
-
-long long
-llroundl (long double x)
-{
-  /* Add +/- 0.5, then round towards zero.  */
-  long double tmp = truncl (x + (x >= 0.0L ?  0.5L : -0.5L));
-  if (!isfinite (tmp) 
-      || tmp > (long double)LONG_LONG_MAX
-      || tmp < (long double)LONG_LONG_MIN)
-    { 
-      errno = ERANGE;
-      /* Undefined behaviour, so we could return anything.  */
-      /* return tmp > 0.0L ? LONG_LONG_MAX : LONG_LONG_MIN; */
-    }
-  return (long long)tmp;
-}
diff --git a/winsup/mingw/mingwex/math/lround.c b/winsup/mingw/mingwex/math/lround.c
deleted file mode 100644
index 7ee50df90..000000000
--- a/winsup/mingw/mingwex/math/lround.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <math.h>
-#include <limits.h>
-#include <errno.h>
-
-long
-lround (double x)
-{
-  /* Add +/- 0.5 then then round towards zero.  */
-  double tmp = trunc (x + (x >= 0.0 ?  0.5 : -0.5));
-  if (!isfinite (tmp) 
-      || tmp > (double)LONG_MAX
-      || tmp < (double)LONG_MIN)
-    { 
-      errno = ERANGE;
-      /* Undefined behaviour, so we could return anything.  */
-      /* return tmp > 0.0 ? LONG_MAX : LONG_MIN;  */
-    }
-  return (long)tmp;
-}
diff --git a/winsup/mingw/mingwex/math/lround_generic.c b/winsup/mingw/mingwex/math/lround_generic.c
new file mode 100644
index 000000000..e0056a665
--- /dev/null
+++ b/winsup/mingw/mingwex/math/lround_generic.c
@@ -0,0 +1,64 @@
+/*
+ * lround_generic.c
+ *
+ * $Id$
+ *
+ * Provides a generic implementation for the `lround()', `lroundf()',
+ * `lroundl()', `llround()', `llroundf()' and `llroundl()' functions;
+ * compile with `-D FUNCTION=name', with `name' set to each of these
+ * six in turn, to create separate object files for each of the six
+ * functions.
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ *
+ * This is free software.  You may redistribute and/or modify it as you
+ * see fit, without restriction of copyright.
+ *
+ * This software is provided "as is", in the hope that it may be useful,
+ * but WITHOUT WARRANTY OF ANY KIND, not even any implied warranty of
+ * MERCHANTABILITY, nor of FITNESS FOR ANY PARTICULAR PURPOSE.  At no
+ * time will the author accept any form of liability for any damages,
+ * however caused, resulting from the use of this software.
+ *
+ */
+#ifndef FUNCTION
+/*
+ * Normally specified with `-D FUNCTION=name', on the command line.
+ * Valid FUNCTION names are `lround', `lroundf', `lroundl', `llround'
+ * `llroundf' and `llroundl'; specifying anything else will most likely
+ * cause a compilation error.  If user did not specify an appropriate
+ * FUNCTION name, default to `lround'.
+ */
+#define FUNCTION lround
+#endif
+
+#include "round_internal.h"
+
+#include <limits.h>
+#include <errno.h>
+
+/* Generic implementation.
+ * The user is required to specify the FUNCTION name;
+ * the RETURN_TYPE and INPUT_TYPE macros resolve to appropriate
+ * type declarations, to match the selected FUNCTION prototype,
+ * while RETURN_MAX and RETURN_MIN map to the correspondingly
+ * appropriate limits.h manifest values, to establish the
+ * valid range for the RETURN_TYPE.
+ */
+RETURN_TYPE FUNCTION( INPUT_TYPE x )
+{
+  if( !isfinite( x ) || !isfinite( x = round_internal( x ) )
+  ||  (x > MAX_RETURN_VALUE) || (x < MIN_RETURN_VALUE)        )
+    /*
+     * Undefined behaviour...
+     * POSIX requires us to report a domain error; ANSI C99 says we
+     * _may_ report a range error, and previous MinGW implementation
+     * set `errno = ERANGE' here; we change that, conforming to the
+     * stricter requiremment of the POSIX standard.
+     */
+    errno = EDOM;
+
+  return (RETURN_TYPE)(x);
+}
+
+/* $RCSfile$$Revision$: end of file */
diff --git a/winsup/mingw/mingwex/math/lroundf.c b/winsup/mingw/mingwex/math/lroundf.c
deleted file mode 100644
index 82df69842..000000000
--- a/winsup/mingw/mingwex/math/lroundf.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <math.h>
-#include <limits.h>
-#include <errno.h>
-
-long
-lroundf (float x)
-{
-  /* Add +/- 0.5, then round towards zero.  */
-  float tmp = truncf (x + (x >= 0.0F ?  0.5F : -0.5F));
-  if (!isfinite (tmp) 
-      || tmp > (float)LONG_MAX
-      || tmp < (float)LONG_MIN)
-    { 
-      errno = ERANGE;
-      /* Undefined behaviour, so we could return anything.  */
-      /* return tmp > 0.0F ? LONG_MAX : LONG_MIN;  */
-    }
-  return (long)tmp;
-}  
diff --git a/winsup/mingw/mingwex/math/lroundl.c b/winsup/mingw/mingwex/math/lroundl.c
deleted file mode 100644
index 7a6348124..000000000
--- a/winsup/mingw/mingwex/math/lroundl.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <math.h>
-#include <limits.h>
-#include <errno.h>
-
-long
-lroundl (long double x)
-{
-  /* Add +/- 0.5, then round towards zero.  */
-  long double tmp = truncl (x + (x >= 0.0L ?  0.5L : -0.5L));
-  if (!isfinite (tmp) 
-      || tmp > (long double)LONG_MAX
-      || tmp < (long double)LONG_MIN)
-    { 
-      errno = ERANGE;
-      /* Undefined behaviour, so we could return anything.  */
-      /* return tmp > 0.0L ? LONG_MAX : LONG_MIN;  */
-    }
-  return (long)tmp;
-}
diff --git a/winsup/mingw/mingwex/math/round.c b/winsup/mingw/mingwex/math/round.c
deleted file mode 100644
index d2d4cab5e..000000000
--- a/winsup/mingw/mingwex/math/round.c
+++ /dev/null
@@ -1,8 +0,0 @@
-#include <math.h>
-
-double
-round (double x)
-{
-  /* Add +/- 0.5 then then round towards zero.  */
-  return trunc ( x + (x >= 0.0 ?  0.5 : -0.5));
-}
diff --git a/winsup/mingw/mingwex/math/round_generic.c b/winsup/mingw/mingwex/math/round_generic.c
new file mode 100644
index 000000000..b3e0f2608
--- /dev/null
+++ b/winsup/mingw/mingwex/math/round_generic.c
@@ -0,0 +1,51 @@
+/*
+ * round_generic.c
+ *
+ * $Id$
+ *
+ * Provides a generic implementation for the `round()', `roundf()'
+ * and `roundl()' functions; compile with `-D FUNCTION=name', with
+ * `name' set to each of these three in turn, to create separate
+ * object files for each of the three functions.
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ *
+ * This is free software.  You may redistribute and/or modify it as you
+ * see fit, without restriction of copyright.
+ *
+ * This software is provided "as is", in the hope that it may be useful,
+ * but WITHOUT WARRANTY OF ANY KIND, not even any implied warranty of
+ * MERCHANTABILITY, nor of FITNESS FOR ANY PARTICULAR PURPOSE.  At no
+ * time will the author accept any form of liability for any damages,
+ * however caused, resulting from the use of this software.
+ *
+ */
+#ifndef FUNCTION
+/*
+ * Normally specified with `-D FUNCTION=name', on the command line.
+ * Valid FUNCTION names are `round', `roundf' and `roundl'; specifying
+ * anything else will most likely cause a compilation error.  If user
+ * did not specify any FUNCTION name, default to `round'.
+ */
+#define FUNCTION round
+#endif
+
+#include "round_internal.h"
+
+/* Generic implementation.
+ * The user is required to specify the FUNCTION name;
+ * the RETURN_TYPE and INPUT_TYPE macros resolve to appropriate
+ * type declarations, to match the selected FUNCTION prototype.
+ */
+RETURN_TYPE FUNCTION( INPUT_TYPE x )
+{
+  /* Round to nearest integer, away from zero for half-way.
+   *
+   * We split it with the `round_internal()' function in
+   * a private header file, so that it may be shared by this,
+   * `lround()' and `llround()' implementations.
+   */
+  return isfinite( x ) ? round_internal( x ) : x;
+}
+
+/* $RCSfile$$Revision$: end of file */
diff --git a/winsup/mingw/mingwex/math/round_internal.h b/winsup/mingw/mingwex/math/round_internal.h
new file mode 100644
index 000000000..bb3c5ad23
--- /dev/null
+++ b/winsup/mingw/mingwex/math/round_internal.h
@@ -0,0 +1,155 @@
+#ifndef _ROUND_INTERNAL_H
+/*
+ * round_internal.h
+ *
+ * $Id$
+ *
+ * Provides a generic implementation of the numerical rounding
+ * algorithm, which is shared by all functions in the `round()',
+ * `lround()' and `llround()' families.
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ *
+ * This is free software.  You may redistribute and/or modify it as you
+ * see fit, without restriction of copyright.
+ *
+ * This software is provided "as is", in the hope that it may be useful,
+ * but WITHOUT WARRANTY OF ANY KIND, not even any implied warranty of
+ * MERCHANTABILITY, nor of FITNESS FOR ANY PARTICULAR PURPOSE.  At no
+ * time will the author accept any form of liability for any damages,
+ * however caused, resulting from the use of this software.
+ *
+ */
+#define _ROUND_INTERNAL_H
+
+#include <math.h>
+#include <fenv.h>
+
+#define TYPE_PASTE( NAME, TYPE )    NAME##TYPE
+
+#define INPUT_TYPE                  INPUT_TYPEDEF( FUNCTION )
+#define INPUT_TYPEDEF( FUNCTION )   TYPE_PASTE( FUNCTION, _input_type )
+/*
+ * The types for the formal parameter, to each of the derived functions.
+ */
+#define round_input_type            double
+#define roundf_input_type           float
+#define roundl_input_type           long double
+
+#define lround_input_type           double
+#define lroundf_input_type          float
+#define lroundl_input_type          long double
+
+#define llround_input_type          double
+#define llroundf_input_type         float
+#define llroundl_input_type         long double
+
+#define RETURN_TYPE                 RETURN_TYPEDEF( FUNCTION )
+#define RETURN_TYPEDEF( FUNCTION )  TYPE_PASTE( FUNCTION, _return_type )
+/*
+ * The types for the return value, from each of the derived functions.
+ */
+#define round_return_type           double
+#define roundf_return_type          float
+#define roundl_return_type          long double
+
+#define lround_return_type          long
+#define lroundf_return_type         long
+#define lroundl_return_type         long
+
+#define llround_return_type         long long
+#define llroundf_return_type        long long
+#define llroundl_return_type        long long
+
+#define MAX_RETURN_VALUE            RETURN_MAX( FUNCTION )
+#define RETURN_MAX( FUNCTION )      TYPE_PASTE( FUNCTION, _return_max )
+/*
+ * The maximum values which may be returned by each of the derived functions
+ * in the `lround' or the `llround' families.
+ */
+#define lround_return_max           LONG_MAX
+#define lroundf_return_max          LONG_MAX
+#define lroundl_return_max          LONG_MAX
+
+#define llround_return_max          LLONG_MAX
+#define llroundf_return_max         LLONG_MAX
+#define llroundl_return_max         LLONG_MAX
+
+#define MIN_RETURN_VALUE            RETURN_MIN( FUNCTION )
+#define RETURN_MIN( FUNCTION )      TYPE_PASTE( FUNCTION, _return_min )
+/*
+ * The minimum values which may be returned by each of the derived functions
+ * in the `lround' or the `llround' families.
+ */
+#define lround_return_min           LONG_MIN
+#define lroundf_return_min          LONG_MIN
+#define lroundl_return_min          LONG_MIN
+
+#define llround_return_min          LLONG_MIN
+#define llroundf_return_min         LLONG_MIN
+#define llroundl_return_min         LLONG_MIN
+
+#define REF_VALUE( VALUE )          REF_TYPE( FUNCTION, VALUE )
+#define REF_TYPE( FUNC, VAL )       TYPE_PASTE( FUNC, _ref )( VAL )
+/*
+ * Macros for expressing constant values of the appropriate data type,
+ * for use in each of the derived functions.
+ */
+#define round_ref( VALUE )          VALUE
+#define lround_ref( VALUE )         VALUE
+#define llround_ref( VALUE )        VALUE
+
+#define roundl_ref( VALUE )         TYPE_PASTE( VALUE, L )
+#define lroundl_ref( VALUE )        TYPE_PASTE( VALUE, L )
+#define llroundl_ref( VALUE )       TYPE_PASTE( VALUE, L )
+
+#define roundf_ref( VALUE )         TYPE_PASTE( VALUE, F )
+#define lroundf_ref( VALUE )        TYPE_PASTE( VALUE, F )
+#define llroundf_ref( VALUE )       TYPE_PASTE( VALUE, F )
+
+static __inline__
+INPUT_TYPE __attribute__(( always_inline )) round_internal( INPUT_TYPE x )
+#define ROUND_MODES ( FE_TONEAREST | FE_UPWARD | FE_DOWNWARD | FE_TOWARDZERO )
+{
+  /* Generic helper function, for rounding of the input parameter value to
+   * the nearest integer value.
+   */
+  INPUT_TYPE z;
+  unsigned short saved_CW, tmp_required_CW;
+
+  /* Rounding method suggested by Danny Smith <dannysmith@users.sf.net>
+   *
+   * Save the FPU control word state, set rounding mode TONEAREST, round the
+   * input value, then restore the original FPU control word state.
+   */
+  __asm__( "fnstcw %0;" : "=m"( saved_CW ));
+  tmp_required_CW = ( saved_CW & ~ROUND_MODES ) | FE_TONEAREST;
+  __asm__( "fldcw %0;" :: "m"( tmp_required_CW ));
+  __asm__( "frndint;" : "=t"( z ) : "0"( x ));
+  __asm__( "fldcw %0;" :: "m"( saved_CW ));
+
+  /* We now have a possible rounded value; unfortunately the FPU uses the
+   * `round-to-even' rule for exact mid-way cases, where both C99 and POSIX
+   * require us to always round away from zero, so we need to adjust those
+   * mid-way cases which the FPU rounded in the wrong direction.
+   *
+   * Correction method suggested by Greg Chicares <gchicares@sbcglobal.net>
+   */
+  return x < REF_VALUE( 0.0 )
+    ? /*
+       * For negative input values, an incorrectly rounded value will be
+       * exactly 0.5 greater than the original value; when we find such an
+       * exact rounding offset, we must subtract an additional 1.0 from the
+       * rounded result, otherwise we return the rounded result unchanged.
+       */
+      z - x == REF_VALUE( 0.5 ) ? z - REF_VALUE( 1.0 ) : z
+
+    : /* For positive input values, an incorrectly rounded value will be
+       * exactly 0.5 less than the original value; when we find such an exact
+       * rounding offset, we must add an additional 1.0 to the rounded result,
+       * otherwise we return the rounded result unchanged.
+       */
+      x - z == REF_VALUE( 0.5 ) ? z + REF_VALUE( 1.0 ) : z;
+}
+
+#endif /* !defined _ROUND_INTERNAL_H: $RCSfile$: end of file */
diff --git a/winsup/mingw/mingwex/math/roundf.c b/winsup/mingw/mingwex/math/roundf.c
deleted file mode 100644
index b50d950a7..000000000
--- a/winsup/mingw/mingwex/math/roundf.c
+++ /dev/null
@@ -1,8 +0,0 @@
-#include <math.h>
-
-float
-roundf (float x)
-{
-  /* Add +/- 0.5 then then round towards zero.  */
-  return truncf ( x + (x >= 0.0F ?  0.5F : -0.5F));
-}
diff --git a/winsup/mingw/mingwex/math/roundl.c b/winsup/mingw/mingwex/math/roundl.c
deleted file mode 100644
index 9c5f0aca1..000000000
--- a/winsup/mingw/mingwex/math/roundl.c
+++ /dev/null
@@ -1,8 +0,0 @@
-#include <math.h>
-
-long double
-roundl (long double x)
-{
-  /* Add +/- 0.5 then then round towards zero.  */
-  return truncl ( x + (x >= 0.0L ?  0.5L : -0.5L));
-}