> 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 <br...@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.
=============================== 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
