master 5e229f8 1/2: Fix rounding errors in time conversion

[Paul Eggert]

branch: master
commit 5e229f88f946c9c246966defeb629965f65d9549
Author: Paul Eggert <address@hidden>
Commit: Paul Eggert <address@hidden>

    Fix rounding errors in time conversion
    
    * src/timefns.c (frac_to_double): Pass FLT_RADIX to mpz_sizeinbase
    instead of doing the radix calculation ourselves, not always
    correctly.  Fix off-by-one error in scale, which caused
    double-rounding.
    (decode_time_components): Use frac_to_double (via decode_ticks_hz)
    to fix double-rounding error that can occur even though
    intermediate results are long double.
    * test/src/timefns-tests.el (float-time-precision):
    Test the above fixes.
---
 src/timefns.c             | 79 +++++++++++++++++++++--------------------------
 test/src/timefns-tests.el |  3 ++
 2 files changed, 38 insertions(+), 44 deletions(-)

diff --git a/src/timefns.c b/src/timefns.c
index b301c58..b96a6cb 100644
--- a/src/timefns.c
+++ b/src/timefns.c
@@ -576,18 +576,14 @@ frac_to_double (Lisp_Object numerator, Lisp_Object 
denominator)
 
   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);
-  ptrdiff_t dbits = mpz_sizeinbase (*d, 2);
-  eassume (0 < nbits);
-  eassume (0 < dbits);
-  ptrdiff_t ndig = (nbits + LOG2_FLT_RADIX - 1) / LOG2_FLT_RADIX;
-  ptrdiff_t ddig = (dbits + LOG2_FLT_RADIX - 1) / LOG2_FLT_RADIX;
+  ptrdiff_t ndig = mpz_sizeinbase (*n, FLT_RADIX);
+  ptrdiff_t ddig = mpz_sizeinbase (*d, FLT_RADIX);
 
   /* Scale with SCALE when doing integer division.  That is, compute
      (N * FLT_RADIX**SCALE) / D [or, if SCALE is negative, N / (D *
      FLT_RADIX**-SCALE)] as a bignum, convert the bignum to double,
      then divide the double by FLT_RADIX**SCALE.  */
-  ptrdiff_t scale = ddig - ndig + DBL_MANT_DIG + 1;
+  ptrdiff_t scale = ddig - ndig + DBL_MANT_DIG;
   if (scale < 0)
     {
       mpz_mul_2exp (mpz[1], *d, - (scale * LOG2_FLT_RADIX));
@@ -615,7 +611,7 @@ frac_to_double (Lisp_Object numerator, Lisp_Object 
denominator)
      round to the nearest integer; otherwise, it is less than
      FLT_RADIX ** (DBL_MANT_DIG + 1) and round it to the nearest
      multiple of FLT_RADIX.  Break ties to even.  */
-  if (mpz_sizeinbase (*q, 2) < DBL_MANT_DIG * LOG2_FLT_RADIX)
+  if (mpz_sizeinbase (*q, FLT_RADIX) <= DBL_MANT_DIG)
     {
       /* Converting to double will use the whole quotient so add 1 to
         its absolute value as per round-to-even; i.e., if the doubled
@@ -746,46 +742,41 @@ decode_time_components (enum timeform form,
   ps = ps % 1000000 + 1000000 * (ps % 1000000 < 0);
   us = us % 1000000 + 1000000 * (us % 1000000 < 0);
 
-  if (result)
+  Lisp_Object hz;
+  switch (form)
     {
-      switch (form)
-       {
-       case TIMEFORM_HI_LO:
-         /* Floats and nil were handled above, so it was an integer.  */
-         mpz_swap (mpz[0], *s);
-         result->hz = make_fixnum (1);
-         break;
-
-       case TIMEFORM_HI_LO_US:
-         mpz_set_ui (mpz[0], us);
-         mpz_addmul_ui (mpz[0], *s, 1000000);
-         result->hz = make_fixnum (1000000);
-         break;
-
-       case TIMEFORM_HI_LO_US_PS:
-         {
-           #if FASTER_TIMEFNS && TRILLION <= ULONG_MAX
-             unsigned long i = us;
-             mpz_set_ui (mpz[0], i * 1000000 + ps);
-             mpz_addmul_ui (mpz[0], *s, TRILLION);
-           #else
-             intmax_t i = us;
-             mpz_set_intmax (mpz[0], i * 1000000 + ps);
-             mpz_addmul (mpz[0], *s, ztrillion);
-           #endif
-           result->hz = trillion;
-         }
-         break;
+    case TIMEFORM_HI_LO:
+      /* Floats and nil were handled above, so it was an integer.  */
+      mpz_swap (mpz[0], *s);
+      hz = make_fixnum (1);
+      break;
 
-       default:
-         eassume (false);
-       }
-      result->ticks = make_integer_mpz ();
+    case TIMEFORM_HI_LO_US:
+      mpz_set_ui (mpz[0], us);
+      mpz_addmul_ui (mpz[0], *s, 1000000);
+      hz = make_fixnum (1000000);
+      break;
+
+    case TIMEFORM_HI_LO_US_PS:
+      {
+       #if FASTER_TIMEFNS && TRILLION <= ULONG_MAX
+         unsigned long i = us;
+         mpz_set_ui (mpz[0], i * 1000000 + ps);
+         mpz_addmul_ui (mpz[0], *s, TRILLION);
+       #else
+         intmax_t i = us;
+         mpz_set_intmax (mpz[0], i * 1000000 + ps);
+         mpz_addmul (mpz[0], *s, ztrillion);
+       #endif
+       hz = trillion;
+      }
+      break;
+
+    default:
+      eassume (false);
     }
-  else
-    *dresult = mpz_get_d (*s) + (us * 1e6L + ps) / 1e12L;
 
-  return 0;
+  return decode_ticks_hz (make_integer_mpz (), hz, result, dresult);
 }
 
 enum { DECODE_SECS_ONLY = WARN_OBSOLETE_TIMESTAMPS + 1 };
diff --git a/test/src/timefns-tests.el b/test/src/timefns-tests.el
index 3967590..b24d443 100644
--- a/test/src/timefns-tests.el
+++ b/test/src/timefns-tests.el
@@ -220,6 +220,9 @@ a fixed place on the right and are padded on the left."
              '(23752 27217))))
 
 (ert-deftest float-time-precision ()
+  (should (= (float-time '(0 1 0 4025)) 1.000000004025))
+  (should (= (float-time '(1000000004025 . 1000000000000)) 1.000000004025))
+
   (should (< 0 (float-time '(1 . 10000000000))))
   (should (< (float-time '(-1 . 10000000000)) 0))