From f48e23eba996566adbe112731ceda4761070fbce Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 3 Mar 2013 04:58:55 -0500 Subject: [PATCH 3/4] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp * libguile/numbers.c (scm_i_big2dbl_2exp): New static function. (scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp', with proper rounding. * test-suite/tests/numbers.test ("inexact->exact"): Add tests. --- libguile/numbers.c | 101 +++++++++++++++-------------------------- test-suite/tests/numbers.test | 45 ++++++++++++++++++ 2 files changed, 81 insertions(+), 65 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index b99a04b..49b4a50 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 8dab29c..e3a1099 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4220,6 +4220,51 @@ (pass-if-exception "nan" exception:out-of-range (inexact->exact +nan.0)) + + (pass-if-equal "round up to odd" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111000101 -> + ;; 11111111111111111111111111111111111111111111111111001000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000101)))) + + (pass-if-equal "round down to odd" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111001011 -> + ;; 11111111111111111111111111111111111111111111111111001000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001011)))) + + (pass-if-equal "round tie up to even" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111011100 -> + ;; 11111111111111111111111111111111111111111111111111100000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b100000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b011100)))) + + (pass-if-equal "round tie down to even" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111000100 -> + ;; 11111111111111111111111111111111111111111111111111000000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000100)))) + + (pass-if-equal "round tie up to next power of two" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111111100 -> + ;; 100000000000000000000000000000000000000000000000000000000 + (expt 2 (+ dbl-mant-dig 3)) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b111100)))) (with-test-prefix "2.0**i to exact and back" (do ((i 0 (1+ i)) -- 1.7.10.4