[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
master bede598 3/3: Fix double-rounding bug in ceiling etc.
From: |
Paul Eggert |
Subject: |
master bede598 3/3: Fix double-rounding bug in ceiling etc. |
Date: |
Wed, 13 Nov 2019 16:10:16 -0500 (EST) |
branch: master
commit bede5984246ba734c93fc28148b5f8e1b14d30c5
Author: Paul Eggert <address@hidden>
Commit: Paul Eggert <address@hidden>
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 7e77dbd..a626845 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 38e1891..1d25add 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 cf634a8..cb04d39 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 62201a8..7f1d469 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)