From 0e1c6ff93f27c939ba9e0df945b16ef98eaaeef1 Mon Sep 17 00:00:00 2001
From: Bruno Haible <bruno@clisp.org>
Date: Sun, 26 Feb 2012 00:36:01 +0100
Subject: [PATCH] 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.
---
 ChangeLog        |  9 +++++++++
 lib/fmodl.c      | 44 ++++++++++++++++++++++++++++++++++++++++++--
 lib/remainder.c  | 23 +++++++++++++++++++++--
 lib/remainderf.c | 23 +++++++++++++++++++++--
 lib/remainderl.c | 23 +++++++++++++++++++++--
 5 files changed, 114 insertions(+), 8 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index e5f0d2001..1b6f2dd5a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,14 @@
 2012-02-25  Bruno Haible  <bruno@clisp.org>
 
+	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.
+
+2012-02-25  Bruno Haible  <bruno@clisp.org>
+
 	Tests for module 'remainderl'.
 	* modules/remainderl-tests: New file.
 	* tests/test-remainderl.c: New file.
diff --git a/lib/fmodl.c b/lib/fmodl.c
index 0c2542fe6..a7ccd8bcd 100644
--- a/lib/fmodl.c
+++ b/lib/fmodl.c
@@ -32,8 +32,48 @@ fmodl (long double x, long double y)
 long double
 fmodl (long double x, long double y)
 {
-  long double i = truncl (x / y);
-  return fmal (- i, y, x);
+  long double q = - truncl (x / y);
+  long double r = fmal (q, y, x); /* = x + q * y, computed in one step */
+  /* Correct possible rounding errors in the quotient x / y.  */
+  if (y >= 0)
+    {
+      if (x >= 0)
+        {
+          /* Expect 0 <= r < y.  */
+          if (r < 0)
+            q += 1, r = fmal (q, y, x);
+          else if (r >= y)
+            q -= 1, r = fmal (q, y, x);
+        }
+      else
+        {
+          /* Expect - y < r <= 0.  */
+          if (r > 0)
+            q -= 1, r = fmal (q, y, x);
+          else if (r <= - y)
+            q += 1, r = fmal (q, y, x);
+        }
+    }
+  else
+    {
+      if (x >= 0)
+        {
+          /* Expect 0 <= r < - y.  */
+          if (r < 0)
+            q -= 1, r = fmal (q, y, x);
+          else if (r >= - y)
+            q += 1, r = fmal (q, y, x);
+        }
+      else
+        {
+          /* Expect y < r <= 0.  */
+          if (r > 0)
+            q += 1, r = fmal (q, y, x);
+          else if (r <= y)
+            q -= 1, r = fmal (q, y, x);
+        }
+    }
+  return r;
 }
 
 #endif
diff --git a/lib/remainder.c b/lib/remainder.c
index c4e7c7ac0..16f09e5db 100644
--- 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;
 }
diff --git a/lib/remainderf.c b/lib/remainderf.c
index 5d9ddb812..d75cf2d39 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
 }
diff --git a/lib/remainderl.c b/lib/remainderl.c
index b43592a64..3476f948e 100644
--- a/lib/remainderl.c
+++ b/lib/remainderl.c
@@ -32,8 +32,27 @@ remainderl (long double x, long double y)
 long double
 remainderl (long double x, long double y)
 {
-  long double i = roundl (x / y);
-  return fmal (- i, y, x);
+  long double q = - roundl (x / y);
+  long double r = fmal (q, y, x); /* = x + q * y, computed in one step */
+  /* Correct possible rounding errors in the quotient x / y.  */
+  long double half_y = 0.5L * y;
+  if (y >= 0)
+    {
+      /* Expect -y/2 <= r <= y/2.  */
+      if (r > half_y)
+        q -= 1, r = fmal (q, y, x);
+      else if (r < - half_y)
+        q += 1, r = fmal (q, y, x);
+    }
+  else
+    {
+      /* Expect y/2 <= r <= -y/2.  */
+      if (r > - half_y)
+        q += 1, r = fmal (q, y, x);
+      else if (r < half_y)
+        q -= 1, r = fmal (q, y, x);
+    }
+  return r;
 }
 
 #endif
-- 
2.11.0