diff --git a/libguile/integers.c b/libguile/integers.c index 102278e44..ebe888d7a 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -416,7 +416,7 @@ scm_integer_floor_remainder_iz (scm_t_inum x, SCM y) SCM scm_integer_floor_remainder_zi (SCM x, scm_t_inum y) { - if (SCM_UNLIKELY (y == 0)) + if (y == 0) scm_num_overflow ("floor-remainder"); else { @@ -1364,3 +1364,142 @@ scm_integer_centered_divide_zz (SCM x, SCM y, SCM *qp, SCM *rp) integer_centered_divide_zz (scm_bignum (x), scm_bignum (y), qp, rp); } +static SCM +integer_round_quotient_zz (struct scm_bignum *x, struct scm_bignum *y) +{ + mpz_t q, r, r2, zx, zy; + int cmp, needs_adjustment; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + mpz_init (q); + mpz_init (r); + mpz_init (r2); + alias_bignum_to_mpz (x, zx); + alias_bignum_to_mpz (y, zy); + + mpz_fdiv_qr (q, r, zx, zy); + mpz_mul_2exp (r2, r, 1); /* r2 = 2*r */ + scm_remember_upto_here_1 (x); + + cmp = mpz_cmpabs (r2, zy); + if (mpz_odd_p (q)) + needs_adjustment = (cmp >= 0); + else + needs_adjustment = (cmp > 0); + scm_remember_upto_here_1 (y); + + if (needs_adjustment) + mpz_add_ui (q, q, 1); + + mpz_clear (r); + mpz_clear (r2); + return take_mpz (q); +} + +SCM +scm_integer_round_quotient_ii (scm_t_inum x, scm_t_inum y) +{ + if (y == 0) + scm_num_overflow ("round-quotient"); + + scm_t_inum q = x / y; + scm_t_inum r = x % y; + scm_t_inum ay = y; + scm_t_inum r2 = 2 * r; + + if (y < 0) + { + ay = -ay; + r2 = -r2; + } + + if (q & 1L) + { + if (r2 >= ay) + q++; + else if (r2 <= -ay) + q--; + } + else + { + if (r2 > ay) + q++; + else if (r2 < -ay) + q--; + } + return long_to_scm (q); +} + +SCM +scm_integer_round_quotient_iz (scm_t_inum x, SCM y) +{ + return integer_round_quotient_zz (long_to_bignum (x), scm_bignum (y)); +} + +SCM +scm_integer_round_quotient_zi (SCM x, scm_t_inum y) +{ + if (y == 0) + scm_num_overflow ("round-quotient"); + if (y == 1) + return x; + + mpz_t q, zx; + mpz_init (q); + alias_bignum_to_mpz (scm_bignum (x), zx); + scm_t_inum r; + int needs_adjustment; + + if (y > 0) + { + r = mpz_fdiv_q_ui (q, zx, y); + if (mpz_odd_p (q)) + needs_adjustment = (2*r >= y); + else + needs_adjustment = (2*r > y); + } + else + { + r = - mpz_cdiv_q_ui (q, zx, -y); + mpz_neg (q, q); + if (mpz_odd_p (q)) + needs_adjustment = (2*r <= y); + else + needs_adjustment = (2*r < y); + } + scm_remember_upto_here_1 (x); + if (needs_adjustment) + mpz_add_ui (q, q, 1); + return take_mpz (q); +} + +SCM +scm_integer_round_quotient_zz (SCM x, SCM y) +{ + SCM q, r, r2; + int cmp, needs_adjustment; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + q = scm_i_mkbig (); + r = scm_i_mkbig (); + r2 = scm_i_mkbig (); + + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ + scm_remember_upto_here_2 (x, r); + + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); + else + needs_adjustment = (cmp > 0); + scm_remember_upto_here_2 (r2, y); + + if (needs_adjustment) + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); + + return scm_i_normbig (q); +} diff --git a/libguile/integers.h b/libguile/integers.h index f941eefc6..dcbc32a10 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -105,6 +105,11 @@ SCM_INTERNAL void scm_integer_centered_divide_zi (SCM x, scm_t_inum y, SCM_INTERNAL void scm_integer_centered_divide_zz (SCM x, SCM y, SCM *qp, SCM *rp); +SCM_INTERNAL SCM scm_integer_round_quotient_ii (scm_t_inum x, scm_t_inum y); +SCM_INTERNAL SCM scm_integer_round_quotient_iz (scm_t_inum x, SCM y); +SCM_INTERNAL SCM scm_integer_round_quotient_zi (SCM x, scm_t_inum y); +SCM_INTERNAL SCM scm_integer_round_quotient_zz (SCM x, SCM y); + #endif /* SCM_INTEGERS_H */ diff --git a/libguile/numbers.c b/libguile/numbers.c index 284d17fd4..d3b794578 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -2439,7 +2439,6 @@ scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp) } static SCM scm_i_inexact_round_quotient (double x, double y); -static SCM scm_i_bigint_round_quotient (SCM x, SCM y); static SCM scm_i_exact_rational_round_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, @@ -2459,55 +2458,15 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_round_quotient { - if (SCM_LIKELY (SCM_I_INUMP (x))) + if (SCM_I_INUMP (x)) { - scm_t_inum xx = SCM_I_INUM (x); - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_round_quotient); - else - { - scm_t_inum qq = xx / yy; - scm_t_inum rr = xx % yy; - scm_t_inum ay = yy; - scm_t_inum r2 = 2 * rr; - - if (SCM_LIKELY (yy < 0)) - { - ay = -ay; - r2 = -r2; - } - - if (qq & 1L) - { - if (r2 >= ay) - qq++; - else if (r2 <= -ay) - qq--; - } - else - { - if (r2 > ay) - qq++; - else if (r2 < -ay) - qq--; - } - if (SCM_LIKELY (SCM_FIXABLE (qq))) - return SCM_I_MAKINUM (qq); - else - return scm_i_inum2big (qq); - } - } + if (SCM_I_INUMP (y)) + return scm_integer_round_quotient_ii (SCM_I_INUM (x), SCM_I_INUM (y)); else if (SCM_BIGP (y)) - { - /* Pass a denormalized bignum version of x (even though it - can fit in a fixnum) to scm_i_bigint_round_quotient */ - return scm_i_bigint_round_quotient (scm_i_long2big (xx), y); - } + return scm_integer_round_quotient_iz (SCM_I_INUM (x), y); else if (SCM_REALP (y)) - return scm_i_inexact_round_quotient (xx, SCM_REAL_VALUE (y)); + return scm_i_inexact_round_quotient (SCM_I_INUM (x), + SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_round_quotient (x, y); else @@ -2516,46 +2475,10 @@ SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, } else if (SCM_BIGP (x)) { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_round_quotient); - else if (SCM_UNLIKELY (yy == 1)) - return x; - else - { - SCM q = scm_i_mkbig (); - scm_t_inum rr; - int needs_adjustment; - - if (yy > 0) - { - rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), yy); - if (mpz_odd_p (SCM_I_BIG_MPZ (q))) - needs_adjustment = (2*rr >= yy); - else - needs_adjustment = (2*rr > yy); - } - else - { - rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), -yy); - mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); - if (mpz_odd_p (SCM_I_BIG_MPZ (q))) - needs_adjustment = (2*rr <= yy); - else - needs_adjustment = (2*rr < yy); - } - scm_remember_upto_here_1 (x); - if (needs_adjustment) - mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); - return scm_i_normbig (q); - } - } + if (SCM_I_INUMP (y)) + return scm_integer_round_quotient_zi (x, SCM_I_INUM (y)); else if (SCM_BIGP (y)) - return scm_i_bigint_round_quotient (x, y); + return scm_integer_round_quotient_zz (x, y); else if (SCM_REALP (y)) return scm_i_inexact_round_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); @@ -2601,38 +2524,6 @@ scm_i_inexact_round_quotient (double x, double y) return scm_i_from_double (scm_c_round (x / y)); } -/* Assumes that both x and y are bigints, though - x might be able to fit into a fixnum. */ -static SCM -scm_i_bigint_round_quotient (SCM x, SCM y) -{ - SCM q, r, r2; - int cmp, needs_adjustment; - - /* Note that x might be small enough to fit into a - fixnum, so we must not let it escape into the wild */ - q = scm_i_mkbig (); - r = scm_i_mkbig (); - r2 = scm_i_mkbig (); - - mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ - scm_remember_upto_here_2 (x, r); - - cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); - if (mpz_odd_p (SCM_I_BIG_MPZ (q))) - needs_adjustment = (cmp >= 0); - else - needs_adjustment = (cmp > 0); - scm_remember_upto_here_2 (r2, y); - - if (needs_adjustment) - mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); - - return scm_i_normbig (q); -} - static SCM scm_i_exact_rational_round_quotient (SCM x, SCM y) {