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.

diff --git a/ChangeLog b/ChangeLog
index 7850616..e7139c0 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,28 @@
 2012-03-04  Bruno Haible  <bruno@clisp.org>
 
+       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.
+
+2012-03-04  Bruno Haible  <bruno@clisp.org>
+
        fmod* tests: More tests.
        * tests/test-fmod.h (my_ldexp): New function.
        (test_function): Reduce amount of random numbers to test. Add tests

diff --git a/lib/remainder.c b/lib/remainder.c
index 16f09e5..869a931 100644 (file)
--- 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;
 }

diff --git a/lib/remainderf.c b/lib/remainderf.c
index d75cf2d..93fd1a3 100644 (file)
--- a/lib/remainderf.c
+++ b/lib/remainderf.c
@@ -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
-}

diff --git a/lib/remainderl.c b/lib/remainderl.c
index 3476f94..474ebe2 100644 (file)
--- a/lib/remainderl.c
+++ b/lib/remainderl.c
@@ -29,30 +29,7 @@ remainderl (long double x, long double y)
 
 #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

diff --git a/m4/remainder.m4 b/m4/remainder.m4
index 855c1bb..f788a20 100644 (file)
--- a/m4/remainder.m4
+++ b/m4/remainder.m4
@@ -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 @@ int main (int argc, char *argv[])
   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])

diff --git a/m4/remainderf.m4 b/m4/remainderf.m4
index 4332bb9..19760ed 100644 (file)
--- a/m4/remainderf.m4
+++ b/m4/remainderf.m4
@@ -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 @@ int main (int argc, char *argv[])
         [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

diff --git a/m4/remainderl.m4 b/m4/remainderl.m4
index 4526bfd..60653bd 100644 (file)
--- a/m4/remainderl.m4
+++ b/m4/remainderl.m4
@@ -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 @@ int main (int argc, char *argv[])
     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

diff --git a/modules/remainder b/modules/remainder
index 342540b..721707f 100644 (file)
--- a/modules/remainder
+++ b/modules/remainder
@@ -8,8 +8,12 @@ m4/mathfunc.m4
 
 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

diff --git a/modules/remainderf b/modules/remainderf
index 69a6426..e4f0600 100644 (file)
--- a/modules/remainderf
+++ b/modules/remainderf
@@ -3,14 +3,19 @@ remainderf() function: floating-point remainder function.
 
 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

diff --git a/modules/remainderl b/modules/remainderl
index 20e2498..4eb1c91 100644 (file)
--- a/modules/remainderl
+++ b/modules/remainderl
@@ -3,14 +3,20 @@ remainderl() function: floating-point remainder function.
 
 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