* 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>
+
+ 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.
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;
}
/* 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
-}
#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
-# 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,
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])
-# 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,
[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
-# 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,
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
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
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
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