From 42ecb421749b7b5ab64943c402c92587509e9f04 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Wed, 23 Feb 2011 00:04:12 -0500 Subject: [PATCH 4/7] Introduce double_precision and use it * libguile/numbers.c (double_precision): New internal static global which contains the number of bits of precision in the C double type. (scm_i_dbl2num, scm_num_eq_p): Use double_precision instead of DBL_MANT_DIG. Previously, the code improperly assumed that DBL_MANT_DIG is a number of bits, when in fact it is a number of base-`FLT_RADIX' digits. (log_of_shifted_double): Use double_precision instead of scm_dblprec[0]. (Note that scm_dblprec[0] is used by the idbl2str, and is either 60 or the number of bits of precision of doubles, whichever is less. This limitation cannot be lifted without increasing the size of the fx_per_radix array, which is already over 16 kilobytes.) --- libguile/numbers.c | 64 ++++++++++++++++++++++++++++++++++++--------------- 1 files changed, 45 insertions(+), 19 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index e2c74d2..1274306 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -168,10 +168,40 @@ scm_from_complex_double (complex double z) #endif /* GUILE_I */ +/* double_precision is the smallest positive integer k + such that 1.0 == 1.0 + 2.0^(-k). -static mpz_t z_negative_one; + For floating-point numbers with significands represented + in base 2, it is the number of bits in the significand + (including the most-significant 1 bit which is sometimes + left implicit and not actually stored). + + When the C "double" type corresponds to the IEEE 754 + binary64 format, double_precision should be 53. +*/ +static long double_precision; + +static void +init_double_precision (void) +{ + double sum; + long i; + + for (i = 1; i < 1000000; i++) + { + sum = 1.0 + ldexp (1.0, -i); + if (sum == 1.0) + { + double_precision = i; + return; + } + } + scm_syserror ("init_double_precision"); +} +static mpz_t z_negative_one; + /* Clear the `mpz_t' embedded in bignum PTR. */ static void finalize_bignum (GC_PTR ptr, GC_PTR data) @@ -304,10 +334,10 @@ scm_i_dbl2num (double u) /* scm_i_big2dbl() rounds to the closest representable double, in accordance with R5RS exact->inexact. - 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. + The approach is to use mpz_get_d to pick out the high + double_precision 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. Bignums exactly half way between representable doubles are rounded to the next higher absolute value (ie. away from zero). This seems like an @@ -321,7 +351,7 @@ scm_i_dbl2num (double u) 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 + high double_precision 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 @@ -339,16 +369,11 @@ scm_i_big2dbl (SCM b) #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) + if (bits > double_precision) { - size_t shift = bits - DBL_MANT_DIG; - mpz_init2 (tmp, DBL_MANT_DIG); + size_t shift = bits - double_precision; + mpz_init2 (tmp, double_precision); mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift); result = ldexp (mpz_get_d (tmp), shift); mpz_clear (tmp); @@ -363,9 +388,9 @@ scm_i_big2dbl (SCM b) result = mpz_get_d (SCM_I_BIG_MPZ (b)); #endif - if (bits > DBL_MANT_DIG) + if (bits > double_precision) { - unsigned long pos = bits - DBL_MANT_DIG - 1; + unsigned long pos = bits - double_precision - 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))) @@ -6332,7 +6357,7 @@ scm_num_eq_p (SCM x, SCM y) double yy = SCM_REAL_VALUE (y); return scm_from_bool ((double) xx == yy - && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1 + && (double_precision >= SCM_I_FIXNUM_BIT-1 || xx == (scm_t_signed_bits) yy)); } else if (SCM_COMPLEXP (y)) @@ -6386,7 +6411,7 @@ scm_num_eq_p (SCM x, SCM y) /* see comments with inum/real above */ scm_t_signed_bits yy = SCM_I_INUM (y); return scm_from_bool (xx == (double) yy - && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1 + && (double_precision >= SCM_I_FIXNUM_BIT-1 || (scm_t_signed_bits) xx == yy)); } else if (SCM_BIGP (y)) @@ -9519,7 +9544,7 @@ log_of_shifted_double (double x, long shift) static SCM log_of_exact_integer_with_size (SCM n, long size) { - long shift = size - 2 * scm_dblprec[0]; + long shift = size - 2 * double_precision; if (shift > 0) return log_of_shifted_double @@ -9806,6 +9831,7 @@ scm_init_numbers () flo_log10e = scm_from_double (M_LOG10E); /* determine floating point precision */ + init_double_precision (); for (i=2; i <= SCM_MAX_DBL_RADIX; ++i) { init_dblprec(&scm_dblprec[i-2],i); -- 1.7.2.5