From 3f5f479b73744d939d15aa1a5736e013ded3ea1d Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Mon, 4 Mar 2013 18:46:33 -0500 Subject: [PATCH 08/10] Improve inexact division of exact integers. * libguile/numbers.c (scm_i_divide2double): New function. (scm_i_divide2double_lo2b): New variable. (scm_i_fraction2double, log_of_fraction): Use 'scm_i_divide2double'. (do_divide): Removed. Its code is now in 'scm_divide'. (scm_divide2real): Removed. Superceded by 'scm_i_divide2double'. (scm_divide): Inherit code from 'do_divide', but without support for forcing a 'double' result (that functionality is now implemented by 'scm_i_divide2double'). Add FIXME comments in cases where divisions might not be as precise as they should be. (scm_init_numbers): Initialize 'scm_i_divide2double_lo2b'. --- libguile/numbers.c | 257 ++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 177 insertions(+), 80 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index 5feed70..601b2a6 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -410,9 +410,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: @@ -468,11 +465,122 @@ scm_i_make_ratio (SCM numerator, SCM denominator) } #undef FUNC_NAME +static mpz_t scm_i_divide2double_lo2b; + +/* Return n/d as a double, where n and d are exact integers. */ +static double +scm_i_divide2double (SCM n, SCM d) +{ + int neg; + mpz_t nn, dd, lo, hi, x; + ssize_t e; + + if (SCM_I_INUMP (d)) + { + if (SCM_UNLIKELY (scm_is_eq (d, SCM_INUM0))) + { + if (scm_is_true (scm_positive_p (n))) + return 1.0 / 0.0; + else if (scm_is_true (scm_negative_p (n))) + return -1.0 / 0.0; + else + return 0.0 / 0.0; + } + mpz_init_set_si (dd, SCM_I_INUM (d)); + } + else + mpz_init_set (dd, SCM_I_BIG_MPZ (d)); + + if (SCM_I_INUMP (n)) + mpz_init_set_si (nn, SCM_I_INUM (n)); + else + mpz_init_set (nn, SCM_I_BIG_MPZ (n)); + + neg = (mpz_sgn (nn) < 0) ^ (mpz_sgn (dd) < 0); + mpz_abs (nn, nn); + mpz_abs (dd, dd); + + /* Now we need to find the value of e such that: + + For e <= 0: + b^{p-1} - 1/2b <= b^-e n / d < b^p - 1/2 + (2 b^p - 1) <= 2 b b^-e n / d < (2 b^p - 1) b + (2 b^p - 1) d <= 2 b b^-e n < (2 b^p - 1) d b + + For e >= 0: + b^{p-1} - 1/2b <= n / b^e d < b^p - 1/2 + (2 b^p - 1) <= 2 b n / b^e d < (2 b^p - 1) b + (2 b^p - 1) d b^e <= 2 b n < (2 b^p - 1) d b b^e + + p == DBL_MANT_DIG + b == FLT_RADIX (assumed 2 here) */ + + /* Initial guess for e */ + e = mpz_sizeinbase (nn, 2) - mpz_sizeinbase (dd, 2) - (DBL_MANT_DIG-1); + if (e < DBL_MIN_EXP - DBL_MANT_DIG) + e = DBL_MIN_EXP - DBL_MANT_DIG; + + /* Compute initial lo, hi, and x */ + mpz_inits (lo, hi, x, NULL); + mpz_mul_2exp (x, nn, 2 + ((e < 0) ? -e : 0)); + mpz_mul (lo, dd, scm_i_divide2double_lo2b); + if (e > 0) + mpz_mul_2exp (lo, lo, e); + mpz_mul_2exp (hi, lo, 1); + + /* Adjust e as needed */ + while (mpz_cmp (x, lo) < 0 && e > DBL_MIN_EXP - DBL_MANT_DIG) + { + mpz_mul_2exp (x, x, 1); + e--; + } + while (mpz_cmp (x, hi) >= 0) + { + mpz_mul_2exp (hi, hi, 1); + e++; + } + + /* Now compute the rounded mantissa: + n / b^e d (if e >= 0) + n b^-e / d (if e <= 0) */ + { + int cmp, needs_adjustment; + double result; + + if (e < 0) + mpz_mul_2exp (nn, nn, -e); + else + mpz_mul_2exp (dd, dd, e); + + /* mpz does not directly support rounded right + shifts, so we have to do it the hard way. + hi == quotient, lo == remainder */ + mpz_fdiv_qr (hi, lo, nn, dd); + mpz_mul_2exp (lo, lo, 1); + + cmp = mpz_cmp (lo, dd); + if (mpz_odd_p (hi)) + needs_adjustment = (cmp >= 0); + else + needs_adjustment = (cmp > 0); + + if (needs_adjustment) + mpz_add_ui (hi, hi, 1); + + if (neg) + mpz_neg (hi, hi); + + result = ldexp (mpz_get_d (hi), e); + mpz_clears (nn, dd, lo, hi, x, NULL); + return result; + } +} + 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 @@ -7965,8 +8073,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; @@ -7985,18 +8093,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); @@ -8046,11 +8146,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; @@ -8061,11 +8157,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); @@ -8074,6 +8166,9 @@ do_divide (SCM x, SCM y, int inexact) scm_num_overflow (s_divide); else #endif + /* FIXME: Precision may be lost here due to: + (1) The cast from 'scm_t_inum' to 'double' + (2) Double rounding */ return scm_from_double ((double) xx / yy); } else if (SCM_COMPLEXP (y)) @@ -8100,7 +8195,7 @@ do_divide (SCM x, SCM y, int inexact) else if (SCM_FRACTIONP (y)) /* a / b/c = ac / b */ return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), - SCM_FRACTION_NUMERATOR (y)); + SCM_FRACTION_NUMERATOR (y)); else SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); } @@ -8144,43 +8239,24 @@ do_divide (SCM x, SCM y, int inexact) return scm_i_normbig (result); } else - { - if (inexact) - return scm_from_double (scm_i_big2dbl (x) / (double) yy); - else return scm_i_make_ratio (x, y); - } + 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 - { - int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - if (divisible_p) - { - SCM result = scm_i_mkbig (); - mpz_divexact (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - return scm_i_normbig (result); - } - else - return scm_i_make_ratio (x, y); - } + int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + if (divisible_p) + { + SCM result = scm_i_mkbig (); + mpz_divexact (SCM_I_BIG_MPZ (result), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (result); + } + else + return scm_i_make_ratio (x, y); } else if (SCM_REALP (y)) { @@ -8190,6 +8266,8 @@ do_divide (SCM x, SCM y, int inexact) scm_num_overflow (s_divide); else #endif + /* FIXME: Precision may be lost here due to: + (1) scm_i_big2dbl (2) Double rounding */ return scm_from_double (scm_i_big2dbl (x) / yy); } else if (SCM_COMPLEXP (y)) @@ -8199,7 +8277,7 @@ do_divide (SCM x, SCM y, int inexact) } else if (SCM_FRACTIONP (y)) return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), - SCM_FRACTION_NUMERATOR (y)); + SCM_FRACTION_NUMERATOR (y)); else SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); } @@ -8214,10 +8292,16 @@ do_divide (SCM x, SCM y, int inexact) scm_num_overflow (s_divide); else #endif + /* FIXME: Precision may be lost here due to: + (1) The cast from 'scm_t_inum' to 'double' + (2) Double rounding */ return scm_from_double (rx / (double) yy); } else if (SCM_BIGP (y)) { + /* FIXME: Precision may be lost here due to: + (1) The conversion from bignum to double + (2) Double rounding */ double dby = mpz_get_d (SCM_I_BIG_MPZ (y)); scm_remember_upto_here_1 (y); return scm_from_double (rx / dby); @@ -8255,12 +8339,18 @@ do_divide (SCM x, SCM y, int inexact) else #endif { + /* FIXME: Precision may be lost here due to: + (1) The conversion from 'scm_t_inum' to double + (2) Double rounding */ double d = yy; return scm_c_make_rectangular (rx / d, ix / d); } } else if (SCM_BIGP (y)) { + /* FIXME: Precision may be lost here due to: + (1) The conversion from bignum to double + (2) Double rounding */ double dby = mpz_get_d (SCM_I_BIG_MPZ (y)); scm_remember_upto_here_1 (y); return scm_c_make_rectangular (rx / dby, ix / dby); @@ -8294,6 +8384,9 @@ do_divide (SCM x, SCM y, int inexact) } else if (SCM_FRACTIONP (y)) { + /* FIXME: Precision may be lost here due to: + (1) The conversion from fraction to double + (2) Double rounding */ double yy = scm_i_fraction2double (y); return scm_c_make_rectangular (rx / yy, ix / yy); } @@ -8311,12 +8404,12 @@ do_divide (SCM x, SCM y, int inexact) else #endif return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x), - scm_product (SCM_FRACTION_DENOMINATOR (x), y)); + scm_product (SCM_FRACTION_DENOMINATOR (x), y)); } else if (SCM_BIGP (y)) { return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x), - scm_product (SCM_FRACTION_DENOMINATOR (x), y)); + scm_product (SCM_FRACTION_DENOMINATOR (x), y)); } else if (SCM_REALP (y)) { @@ -8326,33 +8419,28 @@ do_divide (SCM x, SCM y, int inexact) scm_num_overflow (s_divide); else #endif + /* FIXME: Precision may be lost here due to: + (1) The conversion from fraction to double + (2) Double rounding */ return scm_from_double (scm_i_fraction2double (x) / yy); } else if (SCM_COMPLEXP (y)) { + /* FIXME: Precision may be lost here due to: + (1) The conversion from fraction to double + (2) Double rounding */ a = scm_i_fraction2double (x); goto complex_div; } else if (SCM_FRACTIONP (y)) return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)), - scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))); + scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))); else SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); } 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 @@ -9617,12 +9705,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); } @@ -9890,6 +9977,16 @@ scm_init_numbers () #endif exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2)); + + { + /* Set scm_i_divide2double_lo2b to (2 b^p - 1) */ + mpz_init_set_ui (scm_i_divide2double_lo2b, 1); + mpz_mul_2exp (scm_i_divide2double_lo2b, + scm_i_divide2double_lo2b, + DBL_MANT_DIG + 1); /* 2 b^p */ + mpz_sub_ui (scm_i_divide2double_lo2b, scm_i_divide2double_lo2b, 1); + } + #include "libguile/numbers.x" } -- 1.7.10.4