fmodl, remainder*: Avoid wrong results due to rounding errors.

* lib/fmodl.c (fmodl): Correct the result if it is not within the
expected bounds.
* lib/remainderf.c (remainderf): Likewise.
* lib/remainder.c (remainder): Likewise.
* lib/remainderl.c (remainderl): Likewise.

1 files changed, 21 insertions, 2 deletions

diff --git a/lib/remainderf.c b/lib/remainderf.c
index 5d9ddb8127..d75cf2d39e 100644
--- a/lib/remainderf.c
+++ b/lib/remainderf.c
@@ -25,7 +25,26 @@ remainderf (float x, float y)
 #if HAVE_REMAINDER
   return (float) remainder ((double) x, (double) y);
 #else
-  float i = roundf (x / y);
-  return fmaf (- i, y, x);
+  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;
 #endif
 }