From: Bruno Haible Date: Sun, 4 Mar 2012 22:01:33 +0000 (+0100) Subject: remainder, remainderf, remainderl: Fix computation for large quotients. X-Git-Tag: v0.1~958 X-Git-Url: http://erislabs.net/gitweb/?a=commitdiff_plain;h=9655379852fdbf3a0797035fc30a4060a2b191cf;p=gnulib.git 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 78506169e..e7139c060 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,28 @@ 2012-03-04 Bruno Haible + 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 + 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 16f09e5db..869a93179 100644 --- 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 . */ -#include +#if ! (defined USE_LONG_DOUBLE || defined USE_FLOAT) +# include +#endif /* Specification. */ #include -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 d75cf2d39..93fd1a364 100644 --- a/lib/remainderf.c +++ b/lib/remainderf.c @@ -19,32 +19,17 @@ /* Specification. */ #include +#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 3476f948e..474ebe217 100644 --- 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 855c1bb3f..f788a20fc 100644 --- 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 4332bb946..19760edf7 100644 --- 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 4526bfd95..60653bd20 100644 --- 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 342540bfa..721707f0d 100644 --- 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 69a642621..e4f0600d7 100644 --- 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 20e249840..4eb1c91bb 100644 --- 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