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

[gnulib.git] / lib / remainder.c

diff --git a/lib/remainder.c b/lib/remainder.c
index c4e7c7a..16f09e5 100644 (file)
--- a/lib/remainder.c
+++ b/lib/remainder.c
@@ -22,6 +22,25 @@
 double
 remainder (double x, double y)
 {
-  double i = round (x / y);
-  return fma (- i, y, x);
+  double q = - round (x / y);
+  double r = fma (q, y, x); /* = x + q * y, computed in one step */
+  /* Correct possible rounding errors in the quotient x / y.  */
+  double half_y = 0.5L * y;
+  if (y >= 0)
+    {
+      /* Expect -y/2 <= r <= y/2.  */
+      if (r > half_y)
+        q -= 1, r = fma (q, y, x);
+      else if (r < - half_y)
+        q += 1, r = fma (q, y, x);
+    }
+  else
+    {
+      /* Expect y/2 <= r <= -y/2.  */
+      if (r > - half_y)
+        q += 1, r = fma (q, y, x);
+      else if (r < half_y)
+        q -= 1, r = fma (q, y, x);
+    }
+  return r;
 }