From d754e3f1de92fa2d33edc610f37ff24b58cddbcd Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Fri, 4 Mar 2011 17:40:56 -0500 Subject: [PATCH 7/7] Improve inexact division of exact integers et al * libguile/numbers.c (scm_to_double_2exp): New internal procedure which converts an arbitrary real number to a normalized C double significand and a base-2 exponent, similar to frexp. This is the basis for making several other inexact numerical operations robust when working with large exact integers and rationals. (scm_divide2real): Removed; superceded by scm_i_divide2double. (scm_i_divide2double): Returns the ratio of two arbitrary real number objects as a C double. (scm_i_fraction2double, log_of_fraction): Use scm_i_divide2double instead of scm_divide2real. (scm_divide, do_divide): The core implementation is now in scm_divide. do_divide has been removed. Previously, do_divide was the basis for both scm_divide and scm_divide2real. Remove the `inexact' flag argument, which (when non-zero) requested an inexact answer. The inexact implementation was flawed in several respects, and is now done by scm_i_divide2double. (scm_i_inum_length): New internal procedure that handles the inum case of scm_integer_length and returns the result as a C int. Used by the inum case of scm_to_double_2exp. (scm_integer_length): Moved inum case to scm_i_inum_length. --- libguile/numbers.c | 170 ++++++++++++++++++++++++++++----------------------- 1 files changed, 93 insertions(+), 77 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index c9be26d..b35f66c 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -401,9 +401,6 @@ scm_i_mpz2num (mpz_t b) } } -/* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */ -static SCM scm_divide2real (SCM x, SCM y); - /* scm_i_make_ratio_already_reduced makes a ratio more efficiently, when the following conditions are known in advance: @@ -459,11 +456,74 @@ scm_i_make_ratio (SCM numerator, SCM denominator) } #undef FUNC_NAME +static int scm_i_inum_length (scm_t_inum n); + +/* scm_to_double_2exp() an inexact frexp for SCM real numbers: it + converts x into a normalized significand and exponent such that + x = 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 x is zero, then the returned exponent and significand are both + zero. */ +static double +scm_to_double_2exp (SCM x, long *expon_p) +{ + if (SCM_I_INUMP (x)) + { + int expon = scm_i_inum_length (abs (SCM_I_INUM (x))); + *expon_p = expon; + return ldexp (SCM_I_INUM (x), -expon); + } + else if (SCM_BIGP (x)) + return scm_i_big2dbl_2exp (x, expon_p); + else if (SCM_REALP (x)) + { + int expon; + double signif = frexp (SCM_REAL_VALUE (x), &expon); + *expon_p = expon; + return signif; + } + else if (SCM_FRACTIONP (x)) + { + long n_expon, d_expon; + double n_signif = scm_to_double_2exp (SCM_FRACTION_NUMERATOR (x), + &n_expon); + double d_signif = scm_to_double_2exp (SCM_FRACTION_DENOMINATOR (x), + &d_expon); + long expon = n_expon - d_expon; + double signif = n_signif / d_signif; + if (signif >= 1.0 || signif <= -1.0) + { + signif /= 2.0; + expon++; + } + *expon_p = expon; + return signif; + } + else + scm_wrong_type_arg ("scm_to_double_2exp", SCM_ARG1, x); +} + +/* Return n/d as a double, where n and d are exact integers. */ +static double +scm_i_divide2double (SCM n, SCM d) +{ + if (SCM_LIKELY (SCM_I_INUMP (n) && SCM_I_INUMP (d))) + return ((double) SCM_I_INUM (n)) / ((double) SCM_I_INUM (d)); + else + { + long n_expon, d_expon; + double n_signif = scm_to_double_2exp (n, &n_expon); + double d_signif = scm_to_double_2exp (d, &d_expon); + return ldexp (n_signif / d_signif, n_expon - d_expon); + } +} + double scm_i_fraction2double (SCM z) { - return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z), - SCM_FRACTION_DENOMINATOR (z))); + return scm_i_divide2double (SCM_FRACTION_NUMERATOR (z), + SCM_FRACTION_DENOMINATOR (z)); } static int @@ -5049,6 +5109,22 @@ static const char scm_ilentab[] = { 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4 }; +/* Return the number of bits necessary to represent n */ +static int scm_i_inum_length (scm_t_inum n) +{ + unsigned long c = 0; + unsigned int l = 4; + + if (n < 0) + n = -1 - n; + while (n) + { + c += 4; + l = scm_ilentab [15 & n]; + n >>= 4; + } + return c - 4 + l; +} SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0, (SCM n), @@ -5065,20 +5141,7 @@ SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0, #define FUNC_NAME s_scm_integer_length { if (SCM_I_INUMP (n)) - { - unsigned long c = 0; - unsigned int l = 4; - scm_t_inum nn = SCM_I_INUM (n); - if (nn < 0) - nn = -1 - nn; - while (nn) - { - c += 4; - l = scm_ilentab [15 & nn]; - nn >>= 4; - } - return SCM_I_MAKINUM (c - 4 + l); - } + return SCM_I_MAKINUM (scm_i_inum_length (SCM_I_INUM (n))); else if (SCM_BIGP (n)) { /* mpz_sizeinbase looks at the absolute value of negatives, whereas we @@ -7933,8 +7996,8 @@ SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1, #define s_divide s_scm_i_divide #define g_divide g_scm_i_divide -static SCM -do_divide (SCM x, SCM y, int inexact) +SCM +scm_divide (SCM x, SCM y) #define FUNC_NAME s_divide { double a; @@ -7953,18 +8016,10 @@ do_divide (SCM x, SCM y, int inexact) scm_num_overflow (s_divide); #endif else - { - if (inexact) - return scm_from_double (1.0 / (double) xx); - else return scm_i_make_ratio_already_reduced (SCM_INUM1, x); - } + return scm_i_make_ratio_already_reduced (SCM_INUM1, x); } else if (SCM_BIGP (x)) - { - if (inexact) - return scm_from_double (1.0 / scm_i_big2dbl (x)); - else return scm_i_make_ratio_already_reduced (SCM_INUM1, x); - } + return scm_i_make_ratio_already_reduced (SCM_INUM1, x); else if (SCM_REALP (x)) { double xx = SCM_REAL_VALUE (x); @@ -8014,11 +8069,7 @@ do_divide (SCM x, SCM y, int inexact) #endif } else if (xx % yy != 0) - { - if (inexact) - return scm_from_double ((double) xx / (double) yy); - else return scm_i_make_ratio (x, y); - } + return scm_i_make_ratio (x, y); else { scm_t_inum z = xx / yy; @@ -8029,11 +8080,7 @@ do_divide (SCM x, SCM y, int inexact) } } else if (SCM_BIGP (y)) - { - if (inexact) - return scm_from_double ((double) xx / scm_i_big2dbl (y)); - else return scm_i_make_ratio (x, y); - } + return scm_i_make_ratio (x, y); else if (SCM_REALP (y)) { double yy = SCM_REAL_VALUE (y); @@ -8090,29 +8137,10 @@ do_divide (SCM x, SCM y, int inexact) else if (yy == 1) return x; else - { - if (inexact) - return scm_from_double (scm_i_big2dbl (x) / (double) yy); - else return scm_i_make_ratio (x, y); - } - } - else if (SCM_BIGP (y)) - { - /* big_x / big_y */ - if (inexact) - { - /* It's easily possible for the ratio x/y to fit a double - but one or both x and y be too big to fit a double, - hence the use of mpq_get_d rather than converting and - dividing. */ - mpq_t q; - *mpq_numref(q) = *SCM_I_BIG_MPZ (x); - *mpq_denref(q) = *SCM_I_BIG_MPZ (y); - return scm_from_double (mpq_get_d (q)); - } - else return scm_i_make_ratio (x, y); } + else if (SCM_BIGP (y)) + return scm_i_make_ratio (x, y); else if (SCM_REALP (y)) { double yy = SCM_REAL_VALUE (y); @@ -8273,17 +8301,6 @@ do_divide (SCM x, SCM y, int inexact) else SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide); } - -SCM -scm_divide (SCM x, SCM y) -{ - return do_divide (x, y, 0); -} - -static SCM scm_divide2real (SCM x, SCM y) -{ - return do_divide (x, y, 1); -} #undef FUNC_NAME @@ -9534,12 +9551,11 @@ log_of_fraction (SCM n, SCM d) log_of_exact_integer (d))); else if (scm_is_false (scm_negative_p (n))) return scm_from_double - (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d)))); + (log1p (scm_i_divide2double (scm_difference (n, d), d))); else return scm_c_make_rectangular - (log1p (scm_to_double (scm_divide2real - (scm_difference (scm_abs (n), d), - d))), + (log1p (scm_i_divide2double (scm_difference (scm_abs (n), d), + d)), M_PI); } -- 1.7.2.5