diff --git a/libguile/numbers.c b/libguile/numbers.c index 1f4b9a84d..0abcb2548 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -100,6 +100,13 @@ typedef scm_t_signed_bits scm_t_inum; #define DOUBLE_IS_POSITIVE_INFINITY(x) (isinf(x) && ((x) > 0)) #define DOUBLE_IS_NEGATIVE_INFINITY(x) (isinf(x) && ((x) < 0)) +/* Test an inum to see if it can be converted to a double without loss + of precision. Note that this will sometimes return 0 even when 1 + could have been returned, e.g. for large powers of 2. It is designed + to be a fast check to optimize common cases. */ +#define INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE(n) \ + (SCM_I_FIXNUM_BIT-1 <= DBL_MANT_DIG \ + || ((n) ^ ((n) >> (SCM_I_FIXNUM_BIT-1))) < (1L << DBL_MANT_DIG)) #if ! HAVE_DECL_MPZ_INITS @@ -506,10 +513,10 @@ scm_i_divide2double (SCM n, SCM d) if (SCM_LIKELY (SCM_I_INUMP (d))) { - if (SCM_LIKELY (SCM_I_INUMP (n) - && (SCM_I_FIXNUM_BIT-1 <= DBL_MANT_DIG - || (SCM_I_INUM (n) < (1L << DBL_MANT_DIG) - && SCM_I_INUM (d) < (1L << DBL_MANT_DIG))))) + if (SCM_LIKELY + (SCM_I_INUMP (n) + && INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE (SCM_I_INUM (n)) + && INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE (SCM_I_INUM (d)))) /* If both N and D can be losslessly converted to doubles, then we can rely on IEEE floating point to do proper rounding much faster than we can. */ diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index eca4536a9..3d25e6a8f 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4021,6 +4021,19 @@ (let ((big (ash 1 4096))) (= 1.0 (exact->inexact (/ (1+ big) big))))) + ;; In guile 2.0.9, 'exact->inexact' guaranteed proper rounding when + ;; applied to non-negative fractions, but on 64-bit systems would + ;; sometimes double-round when applied to negative fractions, + ;; specifically when the numerator was a fixnum not exactly + ;; representable as a double. + (with-test-prefix "frac inum/inum, numerator not exactly representable as a double" + (let ((n (+ 1 (expt 2 dbl-mant-dig)))) + (for-each (lambda (d) + (test (/ n d) + (/ n d) + (exact->inexact (/ n d)))) + '(3 5 6 7 9 11 13 17 19 23 0.0 -0.0 +nan.0 +inf.0 -inf.0)))) + (test "round up to odd" ;; ===================================================== ;; 11111111111111111111111111111111111111111111111111000101 ->