1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-05-02 04:40:29 +02:00

Fix rounding in scm_i_divide2double for negative arguments.

* libguile/numbers.c (INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE):
  New macro.
  (scm_i_divide2double): Use INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE to
  determine if our fast path is safe.  Previously, negative arguments
  were not checked properly.

* test-suite/tests/numbers.test (exact->inexact): Add tests.
This commit is contained in:
Mark H Weaver 2013-07-16 00:00:23 -04:00
parent 7e8166f5bd
commit 4cc2e41cf7
2 changed files with 24 additions and 4 deletions

View file

@ -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_POSITIVE_INFINITY(x) (isinf(x) && ((x) > 0))
#define DOUBLE_IS_NEGATIVE_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 #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 (d)))
{ {
if (SCM_LIKELY (SCM_I_INUMP (n) if (SCM_LIKELY
&& (SCM_I_FIXNUM_BIT-1 <= DBL_MANT_DIG (SCM_I_INUMP (n)
|| (SCM_I_INUM (n) < (1L << DBL_MANT_DIG) && INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE (SCM_I_INUM (n))
&& SCM_I_INUM (d) < (1L << DBL_MANT_DIG))))) && INUM_LOSSLESSLY_CONVERTIBLE_TO_DOUBLE (SCM_I_INUM (d))))
/* If both N and D can be losslessly converted to doubles, then /* If both N and D can be losslessly converted to doubles, then
we can rely on IEEE floating point to do proper rounding much we can rely on IEEE floating point to do proper rounding much
faster than we can. */ faster than we can. */

View file

@ -4021,6 +4021,19 @@
(let ((big (ash 1 4096))) (let ((big (ash 1 4096)))
(= 1.0 (exact->inexact (/ (1+ big) big))))) (= 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" (test "round up to odd"
;; ===================================================== ;; =====================================================
;; 11111111111111111111111111111111111111111111111111000101 -> ;; 11111111111111111111111111111111111111111111111111000101 ->