From bede5984246ba734c93fc28148b5f8e1b14d30c5 Mon Sep 17 00:00:00 2001 From: Paul Eggert Date: Wed, 13 Nov 2019 13:07:01 -0800 Subject: [PATCH] Fix double-rounding bug in ceiling etc. This is doable now that we have bignums. * src/floatfns.c (integer_value): Remove; no longer used. (rescale_for_division): New function. (rounding_driver): Use it to divide properly (by using bignums) even when arguments are float, fixing a double-rounding FIXME. * src/lisp.h (LOG2_FLT_RADIX): Move here ... * src/timefns.c (frac_to_double): ... from here. * test/src/floatfns-tests.el (big-round): Add a test to catch the double-rounding bug. --- src/floatfns.c | 115 ++++++++++++++++--------------------- src/lisp.h | 2 + src/timefns.c | 2 - test/src/floatfns-tests.el | 4 +- 4 files changed, 55 insertions(+), 68 deletions(-) diff --git a/src/floatfns.c b/src/floatfns.c index 7e77dbd16dc..a626845377a 100644 --- a/src/floatfns.c +++ b/src/floatfns.c @@ -335,19 +335,6 @@ This is the same as the exponent of a float. */) return make_fixnum (value); } -/* True if A is exactly representable as an integer. */ - -static bool -integer_value (Lisp_Object a) -{ - if (FLOATP (a)) - { - double d = XFLOAT_DATA (a); - return d == floor (d) && isfinite (d); - } - return true; -} - /* Return the integer exponent E such that D * FLT_RADIX**E (i.e., scalbn (D, E)) is an integer that has precision equal to D and is representable as a double. @@ -371,70 +358,68 @@ double_integer_scale (double d) && (FP_ILOGB0 != INT_MIN || d != 0))))); } +/* Convert the Lisp number N to an integer and return a pointer to the + converted integer, represented as an mpz_t *. Use *T as a + temporary; the returned value might be T. Scale N by the maximum + of NSCALE and DSCALE while converting. If NSCALE is nonzero, N + must be a float; signal an overflow if NSCALE is greater than + DBL_MANT_DIG - DBL_MIN_EXP, otherwise scalbn (XFLOAT_DATA (N), NSCALE) + must return an integer value, without rounding or overflow. */ + +static mpz_t const * +rescale_for_division (Lisp_Object n, mpz_t *t, int nscale, int dscale) +{ + mpz_t const *pn; + + if (FLOATP (n)) + { + if (DBL_MANT_DIG - DBL_MIN_EXP < nscale) + overflow_error (); + mpz_set_d (*t, scalbn (XFLOAT_DATA (n), nscale)); + pn = t; + } + else + pn = bignum_integer (t, n); + + if (nscale < dscale) + { + emacs_mpz_mul_2exp (*t, *pn, (dscale - nscale) * LOG2_FLT_RADIX); + pn = t; + } + return pn; +} + /* the rounding functions */ static Lisp_Object -rounding_driver (Lisp_Object arg, Lisp_Object divisor, +rounding_driver (Lisp_Object n, Lisp_Object d, double (*double_round) (double), void (*int_divide) (mpz_t, mpz_t const, mpz_t const), EMACS_INT (*fixnum_divide) (EMACS_INT, EMACS_INT)) { - CHECK_NUMBER (arg); + CHECK_NUMBER (n); - double d; - if (NILP (divisor)) - { - if (! FLOATP (arg)) - return arg; - d = XFLOAT_DATA (arg); - } - else - { - CHECK_NUMBER (divisor); - if (integer_value (arg) && integer_value (divisor)) - { - /* Divide as integers. Converting to double might lose - info, even for fixnums; also see the FIXME below. */ - - if (FLOATP (arg)) - arg = double_to_integer (XFLOAT_DATA (arg)); - if (FLOATP (divisor)) - divisor = double_to_integer (XFLOAT_DATA (divisor)); - - if (FIXNUMP (divisor)) - { - if (XFIXNUM (divisor) == 0) - xsignal0 (Qarith_error); - if (FIXNUMP (arg)) - return make_int (fixnum_divide (XFIXNUM (arg), - XFIXNUM (divisor))); - } - int_divide (mpz[0], - *bignum_integer (&mpz[0], arg), - *bignum_integer (&mpz[1], divisor)); - return make_integer_mpz (); - } + if (NILP (d)) + return FLOATP (n) ? double_to_integer (double_round (XFLOAT_DATA (n))) : n; - double f1 = XFLOATINT (arg); - double f2 = XFLOATINT (divisor); - if (! IEEE_FLOATING_POINT && f2 == 0) - xsignal0 (Qarith_error); - /* FIXME: This division rounds, so the result is double-rounded. */ - d = f1 / f2; - } + CHECK_NUMBER (d); - /* Round, coarsely test for fixnum overflow before converting to - EMACS_INT (to avoid undefined C behavior), and then exactly test - for overflow after converting (as FIXNUM_OVERFLOW_P is inaccurate - on floats). */ - double dr = double_round (d); - if (fabs (dr) < 2 * (MOST_POSITIVE_FIXNUM + 1)) + if (FIXNUMP (d)) { - EMACS_INT ir = dr; - if (! FIXNUM_OVERFLOW_P (ir)) - return make_fixnum (ir); + if (XFIXNUM (d) == 0) + xsignal0 (Qarith_error); + + /* Divide fixnum by fixnum specially, for speed. */ + if (FIXNUMP (n)) + return make_int (fixnum_divide (XFIXNUM (n), XFIXNUM (d))); } - return double_to_integer (dr); + + int nscale = FLOATP (n) ? double_integer_scale (XFLOAT_DATA (n)) : 0; + int dscale = FLOATP (d) ? double_integer_scale (XFLOAT_DATA (d)) : 0; + int_divide (mpz[0], + *rescale_for_division (n, &mpz[0], nscale, dscale), + *rescale_for_division (d, &mpz[1], dscale, nscale)); + return make_integer_mpz (); } static EMACS_INT diff --git a/src/lisp.h b/src/lisp.h index 38e1891c894..1d25add9287 100644 --- a/src/lisp.h +++ b/src/lisp.h @@ -3684,6 +3684,8 @@ extern Lisp_Object string_make_unibyte (Lisp_Object); extern void syms_of_fns (void); /* Defined in floatfns.c. */ +verify (FLT_RADIX == 2 || FLT_RADIX == 16); +enum { LOG2_FLT_RADIX = FLT_RADIX == 2 ? 1 : 4 }; int double_integer_scale (double); #ifndef HAVE_TRUNC extern double trunc (double); diff --git a/src/timefns.c b/src/timefns.c index cf634a82b4e..cb04d39aa13 100644 --- a/src/timefns.c +++ b/src/timefns.c @@ -574,8 +574,6 @@ frac_to_double (Lisp_Object numerator, Lisp_Object denominator) && integer_to_intmax (numerator, &intmax_numerator)) return intmax_numerator; - verify (FLT_RADIX == 2 || FLT_RADIX == 16); - enum { LOG2_FLT_RADIX = FLT_RADIX == 2 ? 1 : 4 }; mpz_t const *n = bignum_integer (&mpz[0], numerator); mpz_t const *d = bignum_integer (&mpz[1], denominator); ptrdiff_t nbits = mpz_sizeinbase (*n, 2); diff --git a/test/src/floatfns-tests.el b/test/src/floatfns-tests.el index 62201a877d0..7f1d4691bf0 100644 --- a/test/src/floatfns-tests.el +++ b/test/src/floatfns-tests.el @@ -122,6 +122,8 @@ (ert-deftest big-round () (should (= (floor 54043195528445955 3) - (floor 54043195528445955 3.0)))) + (floor 54043195528445955 3.0))) + (should (= (floor 1.7976931348623157e+308 5e-324) + (ash (1- (ash 1 53)) 2045)))) (provide 'floatfns-tests) -- 2.39.5