Commit bede5984 authored by Paul Eggert's avatar Paul Eggert

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.
parent 02e637ec
Pipeline #4024 failed with stage
in 61 minutes and 8 seconds
......@@ -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
......
......@@ -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);
......
......@@ -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);
......
......@@ -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)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment