From 5e229f88f946c9c246966defeb629965f65d9549 Mon Sep 17 00:00:00 2001 From: Paul Eggert Date: Tue, 3 Mar 2020 10:17:34 -0800 Subject: [PATCH] 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 b301c58df5a..b96a6cb8875 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 396759023e3..b24d44335e5 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)) -- 2.39.2