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, 79 insertions, 19 deletions

diff --git a/lib/remainder.c b/lib/remainder.c
index 16f09e5db0..869a93179a 100644
--- a/lib/remainder.c
+++ b/lib/remainder.c
@@ -14,33 +14,93 @@
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
 
-#include <config.h>
+#if ! (defined USE_LONG_DOUBLE || defined USE_FLOAT)
+# include <config.h>
+#endif
 
 /* Specification.  */
 #include <math.h>
 
-double
-remainder (double x, double y)
+#ifdef USE_LONG_DOUBLE
+# define REMAINDER remainderl
+# define DOUBLE long double
+# define L_(literal) literal##L
+# define FABS fabsl
+# define FMOD fmodl
+# define ISNAN isnanl
+#elif ! defined USE_FLOAT
+# define REMAINDER remainder
+# define DOUBLE double
+# define L_(literal) literal
+# define FABS fabs
+# define FMOD fmod
+# define ISNAN isnand
+#else /* defined USE_FLOAT */
+# define REMAINDER remainderf
+# define DOUBLE float
+# define L_(literal) literal##f
+# define FABS fabsf
+# define FMOD fmodf
+# define ISNAN isnanf
+#endif
+
+#undef NAN
+#if defined _MSC_VER
+static DOUBLE zero;
+# define NAN (zero / zero)
+#else
+# define NAN (L_(0.0) / L_(0.0))
+#endif
+
+DOUBLE
+REMAINDER (DOUBLE x, DOUBLE y)
 {
-  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)
+  if (isfinite (x) && isfinite (y) && y != L_(0.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);
+      if (x == L_(0.0))
+        /* Return x, regardless of the sign of y.  */
+        return x;
+
+      {
+        int negate = ((!signbit (x)) ^ (!signbit (y)));
+        DOUBLE r;
+
+        /* Take the absolute value of x and y.  */
+        x = FABS (x);
+        y = FABS (y);
+
+        /* Trivial case that requires no computation.  */
+        if (x <= L_(0.5) * y)
+          return (negate ? - x : x);
+
+        /* With a fixed y, the function x -> remainder(x,y) has a period 2*y.
+           Therefore we can reduce the argument x modulo 2*y.  And it's no
+           problem if 2*y overflows, since fmod(x,Inf) = x.  */
+        x = FMOD (x, L_(2.0) * y);
+
+        /* Consider the 3 cases:
+             0 <= x <= 0.5 * y
+             0.5 * y < x < 1.5 * y
+             1.5 * y <= x <= 2.0 * y  */
+        if (x <= L_(0.5) * y)
+          r = x;
+        else
+          {
+            r = x - y;
+            if (r > L_(0.5) * y)
+              r = x - L_(2.0) * y;
+          }
+        return (negate ? - r : r);
+      }
     }
   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);
+      if (ISNAN (x) || ISNAN (y))
+        return x + y; /* NaN */
+      else if (isinf (y))
+        return x;
+      else
+        /* x infinite or y zero */
+        return NAN;
     }
-  return r;
 }