master bede598 3/3: Fix double-rounding bug in ceiling etc.

[Paul Eggert]

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)