remainder, remainderf, remainderl: Fix computation for large quotients.

* lib/remainder.c: Completely rewritten.
* lib/remainderf.c (remainderf): Use implementation of remainder.c with
USE_FLOAT.
* lib/remainderl.c (remainderl): Use implementation of remainder.c with
USE_LONG_DOUBLE.
* modules/remainder (Depends-on): Add isfinite, signbit, fabs, fmod,
isnand, isinf. Remove round, fma.
* modules/remainderf (Files): Add lib/remainder.c.
(Depends-on): Add isfinite, signbit, fabsf, fmodf, isnanf, isinf.
Remove roundf, fmaf.
* modules/remainderl (Files): Add lib/remainder.c.
(Depends-on): Add float, isfinite, signbit, fabsl, fmodl, isnanl,
isinf. Remove roundl, fmal.
* m4/remainder.m4 (gl_FUNC_REMAINDER): Update computation of
REMAINDER_LIBM.
* m4/remainderf.m4 (gl_FUNC_REMAINDERF): Update computation of
REMAINDERF_LIBM.
* m4/remainderl.m4 (gl_FUNC_REMAINDERL): Update computation of
REMAINDERL_LIBM.

1 files changed, 8 insertions, 23 deletions

diff --git a/lib/remainderf.c b/lib/remainderf.c
index d75cf2d39e..93fd1a3645 100644
--- a/lib/remainderf.c
+++ b/lib/remainderf.c
@@ -19,32 +19,17 @@
 /* Specification.  */
 #include <math.h>
 
+#if HAVE_REMAINDER
+
 float
 remainderf (float x, float y)
 {
-#if HAVE_REMAINDER
   return (float) remainder ((double) x, (double) y);
+}
+
 #else
-  float q = - roundf (x / y);
-  float r = fmaf (q, y, x); /* = x + q * y, computed in one step */
-  /* Correct possible rounding errors in the quotient x / y.  */
-  float half_y = 0.5L * y;
-  if (y >= 0)
-    {
-      /* Expect -y/2 <= r <= y/2.  */
-      if (r > half_y)
-        q -= 1, r = fmaf (q, y, x);
-      else if (r < - half_y)
-        q += 1, r = fmaf (q, y, x);
-    }
-  else
-    {
-      /* Expect y/2 <= r <= -y/2.  */
-      if (r > - half_y)
-        q += 1, r = fmaf (q, y, x);
-      else if (r < half_y)
-        q -= 1, r = fmaf (q, y, x);
-    }
-  return r;
+
+# define USE_FLOAT
+# include "remainder.c"
+
 #endif
-}