Mercurial > hg > octave-lojdl > gnulib-hg

changeset 16568:8fec0d45d1a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Bruno Haible <bruno@clisp.org>
date Sun, 04 Mar 2012 23:01:33 +0100
parents 05507ae4f78f
children 1f283a634cf1
files ChangeLog lib/remainder.c lib/remainderf.c lib/remainderl.c m4/remainder.m4 m4/remainderf.m4 m4/remainderl.m4 modules/remainder modules/remainderf modules/remainderl
diffstat 10 files changed, 178 insertions(+), 100 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,26 @@
+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.
--- 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;
 }
--- 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
-}
--- a/lib/remainderl.c
+++ b/lib/remainderl.c
@@ -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
--- 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 @@
   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])
--- 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 @@
         [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
--- 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 @@
     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
--- a/modules/remainder
+++ b/modules/remainder
@@ -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
--- a/modules/remainderf
+++ b/modules/remainderf
@@ -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
--- a/modules/remainderl
+++ b/modules/remainderl
@@ -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