Mercurial > hg > octave-lyh > gnulib-hg
changeset 16568:8fec0d45d1a8
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