remainder, remainderf, remainderl: rewrite

[Bruno Haible]

> The additional unit tests for fmod() and fmodl() revealed that gnulib's
> replacement algorithm returned wrong results when the quotient of the
> two arguments was large. For example:
>     fmod (1e30, 1.5)

The gnulib implementation of the remainder() function family had the same
bug. Corrected as follows.


2012-03-04  Bruno Haible  <address@hidden>

        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.

=============================== lib/remainder.c ===============================
/* Remainder.
   Copyright (C) 2012 Free Software Foundation, Inc.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.  */

#if ! (defined USE_LONG_DOUBLE || defined USE_FLOAT)
# include <config.h>
#endif

/* Specification.  */
#include <math.h>

#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)
{
  if (isfinite (x) && isfinite (y) && y != L_(0.0))
    {
      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
    {
      if (ISNAN (x) || ISNAN (y))
        return x + y; /* NaN */
      else if (isinf (y))
        return x;
      else
        /* x infinite or y zero */
        return NAN;
    }
}
===============================================================================
--- lib/remainderf.c.orig       Sun Mar  4 22:56:28 2012
+++ lib/remainderf.c    Sun Mar  4 22:05:15 2012
@@ -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
-}
--- lib/remainderl.c.orig       Sun Mar  4 22:56:28 2012
+++ lib/remainderl.c    Sun Mar  4 22:05:15 2012
@@ -29,30 +29,7 @@
 
 #else
 
-long double
-remainderl (long double x, long double y)
-{
-  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;
-}
+# define USE_LONG_DOUBLE
+# include "remainder.c"
 
 #endif
--- m4/remainder.m4.orig        Sun Mar  4 22:56:28 2012
+++ m4/remainder.m4     Sun Mar  4 22:22:05 2012
@@ -1,4 +1,4 @@
-# remainder.m4 serial 2
+# remainder.m4 serial 3
 dnl Copyright (C) 2012 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -104,18 +104,24 @@
   fi
   if test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1; then
     dnl Find libraries needed to link lib/remainder.c.
-    AC_REQUIRE([gl_FUNC_ROUND])
-    AC_REQUIRE([gl_FUNC_FMA])
+    AC_REQUIRE([gl_FUNC_FABS])
+    AC_REQUIRE([gl_FUNC_FMOD])
+    AC_REQUIRE([gl_FUNC_ISNAND])
     REMAINDER_LIBM=
-    dnl Append $ROUND_LIBM to REMAINDER_LIBM, avoiding gratuitous duplicates.
+    dnl Append $FABS_LIBM to REMAINDER_LIBM, avoiding gratuitous duplicates.
     case " $REMAINDER_LIBM " in
-      *" $ROUND_LIBM "*) ;;
-      *) REMAINDER_LIBM="$REMAINDER_LIBM $ROUND_LIBM" ;;
+      *" $FABS_LIBM "*) ;;
+      *) REMAINDER_LIBM="$REMAINDER_LIBM $FABS_LIBM" ;;
     esac
-    dnl Append $FMA_LIBM to REMAINDER_LIBM, avoiding gratuitous duplicates.
+    dnl Append $FMOD_LIBM to REMAINDER_LIBM, avoiding gratuitous duplicates.
     case " $REMAINDER_LIBM " in
-      *" $FMA_LIBM "*) ;;
-      *) REMAINDER_LIBM="$REMAINDER_LIBM $FMA_LIBM" ;;
+      *" $FMOD_LIBM "*) ;;
+      *) REMAINDER_LIBM="$REMAINDER_LIBM $FMOD_LIBM" ;;
+    esac
+    dnl Append $ISNAND_LIBM to REMAINDER_LIBM, avoiding gratuitous duplicates.
+    case " $REMAINDER_LIBM " in
+      *" $ISNAND_LIBM "*) ;;
+      *) REMAINDER_LIBM="$REMAINDER_LIBM $ISNAND_LIBM" ;;
     esac
   fi
   AC_SUBST([REMAINDER_LIBM])
--- m4/remainderf.m4.orig       Sun Mar  4 22:56:28 2012
+++ m4/remainderf.m4    Sun Mar  4 22:24:23 2012
@@ -1,4 +1,4 @@
-# remainderf.m4 serial 2
+# remainderf.m4 serial 3
 dnl Copyright (C) 2012 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -91,18 +91,24 @@
         [Define to 1 if the remainder() function is available in libc or 
libm.])
       REMAINDERF_LIBM="$REMAINDER_LIBM"
     else
-      AC_REQUIRE([gl_FUNC_ROUNDF])
-      AC_REQUIRE([gl_FUNC_FMAF])
+      AC_REQUIRE([gl_FUNC_FABSF])
+      AC_REQUIRE([gl_FUNC_FMODF])
+      AC_REQUIRE([gl_FUNC_ISNANF])
       REMAINDERF_LIBM=
-      dnl Append $ROUNDF_LIBM to REMAINDERF_LIBM, avoiding gratuitous 
duplicates.
+      dnl Append $FABSF_LIBM to REMAINDERF_LIBM, avoiding gratuitous 
duplicates.
       case " $REMAINDERF_LIBM " in
-        *" $ROUNDF_LIBM "*) ;;
-        *) REMAINDERF_LIBM="$REMAINDERF_LIBM $ROUNDF_LIBM" ;;
+        *" $FABSF_LIBM "*) ;;
+        *) REMAINDERF_LIBM="$REMAINDERF_LIBM $FABSF_LIBM" ;;
       esac
-      dnl Append $FMAF_LIBM to REMAINDERF_LIBM, avoiding gratuitous duplicates.
+      dnl Append $FMODF_LIBM to REMAINDERF_LIBM, avoiding gratuitous 
duplicates.
       case " $REMAINDERF_LIBM " in
-        *" $FMAF_LIBM "*) ;;
-        *) REMAINDERF_LIBM="$REMAINDERF_LIBM $FMAF_LIBM" ;;
+        *" $FMODF_LIBM "*) ;;
+        *) REMAINDERF_LIBM="$REMAINDERF_LIBM $FMODF_LIBM" ;;
+      esac
+      dnl Append $ISNANF_LIBM to REMAINDERF_LIBM, avoiding gratuitous 
duplicates.
+      case " $REMAINDERF_LIBM " in
+        *" $ISNANF_LIBM "*) ;;
+        *) REMAINDERF_LIBM="$REMAINDERF_LIBM $ISNANF_LIBM" ;;
       esac
     fi
   fi
--- m4/remainderl.m4.orig       Sun Mar  4 22:56:29 2012
+++ m4/remainderl.m4    Sun Mar  4 22:25:58 2012
@@ -1,4 +1,4 @@
-# remainderl.m4 serial 2
+# remainderl.m4 serial 3
 dnl Copyright (C) 2012 Free Software Foundation, Inc.
 dnl This file is free software; the Free Software Foundation
 dnl gives unlimited permission to copy and/or distribute it,
@@ -88,18 +88,24 @@
     if test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 1; then
       REMAINDERL_LIBM="$REMAINDER_LIBM"
     else
-      AC_REQUIRE([gl_FUNC_ROUNDL])
-      AC_REQUIRE([gl_FUNC_FMAL])
+      AC_REQUIRE([gl_FUNC_FABSL])
+      AC_REQUIRE([gl_FUNC_FMODL])
+      AC_REQUIRE([gl_FUNC_ISNANL])
       REMAINDERL_LIBM=
-      dnl Append $ROUNDL_LIBM to REMAINDERL_LIBM, avoiding gratuitous 
duplicates.
+      dnl Append $FABSL_LIBM to REMAINDERL_LIBM, avoiding gratuitous 
duplicates.
       case " $REMAINDERL_LIBM " in
-        *" $ROUNDL_LIBM "*) ;;
-        *) REMAINDERL_LIBM="$REMAINDERL_LIBM $ROUNDL_LIBM" ;;
+        *" $FABSL_LIBM "*) ;;
+        *) REMAINDERL_LIBM="$REMAINDERL_LIBM $FABSL_LIBM" ;;
       esac
-      dnl Append $FMAL_LIBM to REMAINDERL_LIBM, avoiding gratuitous duplicates.
+      dnl Append $FMODL_LIBM to REMAINDERL_LIBM, avoiding gratuitous 
duplicates.
       case " $REMAINDERL_LIBM " in
-        *" $FMAL_LIBM "*) ;;
-        *) REMAINDERL_LIBM="$REMAINDERL_LIBM $FMAL_LIBM" ;;
+        *" $FMODL_LIBM "*) ;;
+        *) REMAINDERL_LIBM="$REMAINDERL_LIBM $FMODL_LIBM" ;;
+      esac
+      dnl Append $ISNANL_LIBM to REMAINDERL_LIBM, avoiding gratuitous 
duplicates.
+      case " $REMAINDERL_LIBM " in
+        *" $ISNANL_LIBM "*) ;;
+        *) REMAINDERL_LIBM="$REMAINDERL_LIBM $ISNANL_LIBM" ;;
       esac
     fi
   fi
--- modules/remainder.orig      Sun Mar  4 22:56:29 2012
+++ modules/remainder   Sun Mar  4 22:18:34 2012
@@ -8,8 +8,12 @@
 
 Depends-on:
 math
-round           [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
-fma             [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
+isfinite        [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
+signbit         [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
+fabs            [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
+fmod            [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
+isnand          [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
+isinf           [test $HAVE_REMAINDER = 0 || test $REPLACE_REMAINDER = 1]
 
 configure.ac:
 gl_FUNC_REMAINDER
--- modules/remainderf.orig     Sun Mar  4 22:56:29 2012
+++ modules/remainderf  Sun Mar  4 22:55:25 2012
@@ -3,14 +3,19 @@
 
 Files:
 lib/remainderf.c
+lib/remainder.c
 m4/remainderf.m4
 m4/mathfunc.m4
 
 Depends-on:
 math
 remainder       [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
-roundf          [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
-fmaf            [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
+isfinite        [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
+signbit         [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
+fabsf           [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
+fmodf           [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
+isnanf          [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
+isinf           [test $HAVE_REMAINDERF = 0 || test $REPLACE_REMAINDERF = 1]
 
 configure.ac:
 gl_FUNC_REMAINDERF
--- modules/remainderl.orig     Sun Mar  4 22:56:29 2012
+++ modules/remainderl  Sun Mar  4 22:55:33 2012
@@ -3,14 +3,20 @@
 
 Files:
 lib/remainderl.c
+lib/remainder.c
 m4/remainderl.m4
 m4/mathfunc.m4
 
 Depends-on:
 math
 remainder       [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 1]
-roundl          [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
-fmal            [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+float           [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+isfinite        [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+signbit         [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+fabsl           [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+fmodl           [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+isnanl          [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
+isinf           [{ test $HAVE_REMAINDERL = 0 || test $REPLACE_REMAINDERL = 1; 
} && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0]
 
 configure.ac:
 gl_FUNC_REMAINDERL