diff options
| -rw-r--r-- | libguile/numbers.c | 101 | ||||
| -rw-r--r-- | test-suite/tests/numbers.test | 57 |
2 files changed, 80 insertions, 78 deletions
diff --git a/libguile/numbers.c b/libguile/numbers.c index 3f2afdebb..81461792d 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -330,81 +330,52 @@ scm_i_dbl2num (double u) return scm_i_dbl2big (u); } -/* scm_i_big2dbl() rounds to the closest representable double, in accordance - with R5RS exact->inexact. +static SCM round_right_shift_exact_integer (SCM n, long count); - The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits - (ie. truncate towards zero), then adjust to get the closest double by - examining the next lower bit and adding 1 (to the absolute value) if - necessary. +/* scm_i_big2dbl_2exp() is like frexp for bignums: it converts the + bignum b into a normalized significand and exponent such that + b = significand * 2^exponent and 1/2 <= abs(significand) < 1. + The return value is the significand rounded to the closest + representable double, and the exponent is placed into *expon_p. + If b is zero, then the returned exponent and significand are both + zero. */ - Bignums exactly half way between representable doubles are rounded to the - next higher absolute value (ie. away from zero). This seems like an - adequate interpretation of R5RS "numerically closest", and it's easier - and faster than a full "nearest-even" style. - - The bit test must be done on the absolute value of the mpz_t, which means - we need to use mpz_getlimbn. mpz_tstbit is not right, it treats - negatives as twos complement. - - In GMP before 4.2, mpz_get_d rounding was unspecified. It ended up - following the hardware rounding mode, but applied to the absolute - value of the mpz_t operand. This is not what we want so we put the - high DBL_MANT_DIG bits into a temporary. Starting with GMP 4.2 - (released in March 2006) mpz_get_d now always truncates towards zero. - - ENHANCE-ME: The temporary init+clear to force the rounding in GMP - before 4.2 is a slowdown. It'd be faster to pick out the relevant - high bits with mpz_getlimbn. */ - -double -scm_i_big2dbl (SCM b) +static double +scm_i_big2dbl_2exp (SCM b, long *expon_p) { - double result; - size_t bits; - - bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2); - -#if 1 - { - /* For GMP earlier than 4.2, force truncation towards zero */ - - /* FIXME: DBL_MANT_DIG is the number of base-`FLT_RADIX' digits, - _not_ the number of bits, so this code will break badly on a - system with non-binary doubles. */ - - mpz_t tmp; - if (bits > DBL_MANT_DIG) - { - size_t shift = bits - DBL_MANT_DIG; - mpz_init2 (tmp, DBL_MANT_DIG); - mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift); - result = ldexp (mpz_get_d (tmp), shift); - mpz_clear (tmp); - } - else - { - result = mpz_get_d (SCM_I_BIG_MPZ (b)); - } - } -#else - /* GMP 4.2 or later */ - result = mpz_get_d (SCM_I_BIG_MPZ (b)); -#endif + size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2); + size_t shift = 0; if (bits > DBL_MANT_DIG) { - unsigned long pos = bits - DBL_MANT_DIG - 1; - /* test bit number "pos" in absolute value */ - if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS) - & ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS))) + shift = bits - DBL_MANT_DIG; + b = round_right_shift_exact_integer (b, shift); + if (SCM_I_INUMP (b)) { - result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1); + int expon; + double signif = frexp (SCM_I_INUM (b), &expon); + *expon_p = expon + shift; + return signif; } } - scm_remember_upto_here_1 (b); - return result; + { + long expon; + double signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b)); + scm_remember_upto_here_1 (b); + *expon_p = expon + shift; + return signif; + } +} + +/* scm_i_big2dbl() rounds to the closest representable double, + in accordance with R5RS exact->inexact. */ +double +scm_i_big2dbl (SCM b) +{ + long expon; + double signif = scm_i_big2dbl_2exp (b, &expon); + return ldexp (signif, expon); } SCM diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index bb1424853..6b4e08c3a 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -3858,21 +3858,17 @@ ;;; (with-test-prefix "exact->inexact" - + + ;; Test "(exact->inexact n)", expect "want". + (define (test name n want) + (with-test-prefix name + (pass-if-equal "pos" want (exact->inexact n)) + (pass-if-equal "neg" (- want) (exact->inexact (- n))))) + ;; Test "(exact->inexact n)", expect "want". ;; "i" is a index, for diagnostic purposes. (define (try-i i n want) - (with-test-prefix (list i n want) - (with-test-prefix "pos" - (let ((got (exact->inexact n))) - (pass-if "inexact?" (inexact? got)) - (pass-if (list "=" got) (= want got)))) - (set! n (- n)) - (set! want (- want)) - (with-test-prefix "neg" - (let ((got (exact->inexact n))) - (pass-if "inexact?" (inexact? got)) - (pass-if (list "=" got) (= want got)))))) + (test (list i n want) n want)) (with-test-prefix "2^i, no round" (do ((i 0 (1+ i)) @@ -3945,7 +3941,42 @@ ;; convert the num and den to doubles, resulting in infs. (pass-if "frac big/big, exceeding double" (let ((big (ash 1 4096))) - (= 1.0 (exact->inexact (/ (1+ big) big)))))) + (= 1.0 (exact->inexact (/ (1+ big) big))))) + + (test "round up to odd" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111000101 -> + ;; 11111111111111111111111111111111111111111111111111001000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000101) + (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b001000)) + + (test "round down to odd" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111001011 -> + ;; 11111111111111111111111111111111111111111111111111001000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001011) + (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b001000)) + + (test "round tie up to even" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111011100 -> + ;; 11111111111111111111111111111111111111111111111111100000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b011100) + (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b100000)) + + (test "round tie down to even" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111000100 -> + ;; 11111111111111111111111111111111111111111111111111000000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000100) + (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b000000)) + + (test "round tie up to next power of two" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111111100 -> + ;; 100000000000000000000000000000000000000000000000000000000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b111100) + (expt 2.0 (+ dbl-mant-dig 3)))) ;;; ;;; expt |