diff --git a/libguile/integers.c b/libguile/integers.c index 402b2e3c9..e9c3a067c 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -160,28 +160,26 @@ negate_bignum (struct scm_bignum *z) return z; } -static SCM +static struct scm_bignum * make_bignum_1 (int is_negative, mp_limb_t limb) { struct scm_bignum *z = allocate_bignum (1); z->limbs[0] = limb; - return SCM_PACK (is_negative ? negate_bignum(z) : z); + return is_negative ? negate_bignum(z) : z; } -static SCM +static struct scm_bignum * ulong_to_bignum (unsigned long u) { - ASSERT (!SCM_POSFIXABLE (u)); return make_bignum_1 (0, u); }; -static SCM +static struct scm_bignum * long_to_bignum (long i) { if (i > 0) return ulong_to_bignum (i); - ASSERT (!SCM_NEGFIXABLE (i)); return make_bignum_1 (1, long_magnitude (i)); }; @@ -190,7 +188,7 @@ long_to_scm (long i) { if (SCM_FIXABLE (i)) return SCM_I_MAKINUM (i); - return long_to_bignum (i); + return SCM_PACK (long_to_bignum (i)); } static SCM @@ -198,7 +196,7 @@ ulong_to_scm (unsigned long i) { if (SCM_POSFIXABLE (i)) return SCM_I_MAKINUM (i); - return ulong_to_bignum (i); + return SCM_PACK (ulong_to_bignum (i)); } static struct scm_bignum * @@ -998,3 +996,126 @@ scm_integer_truncate_divide_zz (SCM x, SCM y, SCM *qp, SCM *rp) *qp = take_mpz (q); *rp = take_mpz (r); } + +static SCM +integer_centered_quotient_zz (struct scm_bignum *x, struct scm_bignum *y) +{ + mpz_t q, r, min_r, zx, zy; + mpz_init (q); + mpz_init (r); + mpz_init (min_r); + alias_bignum_to_mpz (x, zx); + alias_bignum_to_mpz (y, zy); + + /* Note that x might be small enough to fit into a fixnum, so we must + not let it escape into the wild. */ + + /* min_r will eventually become -abs(y)/2 */ + mpz_tdiv_q_2exp (min_r, zy, 1); + + /* Arrange for rr to initially be non-positive, because that + simplifies the test to see if it is within the needed bounds. */ + if (mpz_sgn (zy) > 0) + { + mpz_cdiv_qr (q, r, zx, zy); + scm_remember_upto_here_2 (x, y); + mpz_neg (min_r, min_r); + if (mpz_cmp (r, min_r) < 0) + mpz_sub_ui (q, q, 1); + } + else + { + mpz_fdiv_qr (q, r, zx, zy); + scm_remember_upto_here_2 (x, y); + if (mpz_cmp (r, min_r) < 0) + mpz_add_ui (q, q, 1); + } + mpz_clear (r); + mpz_clear (min_r); + return take_mpz (q); +} + +SCM +scm_integer_centered_quotient_ii (scm_t_inum x, scm_t_inum y) +{ + if (y == 0) + scm_num_overflow ("centered-quotient"); + + scm_t_inum q = x / y; + scm_t_inum r = x % y; + if (x > 0) + { + if (y > 0) + { + if (r >= (y + 1) / 2) + q++; + } + else + { + if (r >= (1 - y) / 2) + q--; + } + } + else + { + if (y > 0) + { + if (r < -y / 2) + q--; + } + else + { + if (r < y / 2) + q++; + } + } + return long_to_scm (q); +} + +SCM +scm_integer_centered_quotient_iz (scm_t_inum x, SCM y) +{ + return integer_centered_quotient_zz (long_to_bignum (x), + scm_bignum (y)); +} + +SCM +scm_integer_centered_quotient_zi (SCM x, scm_t_inum y) +{ + if (y == 0) + scm_num_overflow ("centered-quotient"); + else if (y == 1) + return x; + else + { + mpz_t q, zx; + mpz_init (q); + alias_bignum_to_mpz (scm_bignum (x), zx); + scm_t_inum r; + /* Arrange for r to initially be non-positive, because that + simplifies the test to see if it is within the needed + bounds. */ + if (y > 0) + { + r = - mpz_cdiv_q_ui (q, zx, y); + scm_remember_upto_here_1 (x); + if (r < -y / 2) + mpz_sub_ui (q, q, 1); + } + else + { + r = - mpz_cdiv_q_ui (q, zx, -y); + scm_remember_upto_here_1 (x); + mpz_neg (q, q); + if (r < y / 2) + mpz_add_ui (q, q, 1); + } + return take_mpz (q); + } +} + +SCM +scm_integer_centered_quotient_zz (SCM x, SCM y) +{ + return integer_centered_quotient_zz (scm_bignum (x), scm_bignum (y)); +} diff --git a/libguile/integers.h b/libguile/integers.h index bd25a4b48..44dc0e376 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -86,6 +86,11 @@ SCM_INTERNAL void scm_integer_truncate_divide_zi (SCM x, scm_t_inum y, SCM_INTERNAL void scm_integer_truncate_divide_zz (SCM x, SCM y, SCM *qp, SCM *rp); +SCM_INTERNAL SCM scm_integer_centered_quotient_ii (scm_t_inum x, scm_t_inum y); +SCM_INTERNAL SCM scm_integer_centered_quotient_iz (scm_t_inum x, SCM y); +SCM_INTERNAL SCM scm_integer_centered_quotient_zi (SCM x, scm_t_inum y); +SCM_INTERNAL SCM scm_integer_centered_quotient_zz (SCM x, SCM y); + #endif /* SCM_INTEGERS_H */ diff --git a/libguile/numbers.c b/libguile/numbers.c index 00491e171..fd3a8c6c7 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -2111,7 +2111,6 @@ scm_i_exact_rational_truncate_divide (SCM x, SCM y, SCM *qp, SCM *rp) } static SCM scm_i_inexact_centered_quotient (double x, double y); -static SCM scm_i_bigint_centered_quotient (SCM x, SCM y); static SCM scm_i_exact_rational_centered_quotient (SCM x, SCM y); SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, @@ -2129,58 +2128,16 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_centered_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_centered_quotient); - else - { - scm_t_inum qq = xx / yy; - scm_t_inum rr = xx % yy; - if (SCM_LIKELY (xx > 0)) - { - if (SCM_LIKELY (yy > 0)) - { - if (rr >= (yy + 1) / 2) - qq++; - } - else - { - if (rr >= (1 - yy) / 2) - qq--; - } - } - else - { - if (SCM_LIKELY (yy > 0)) - { - if (rr < -yy / 2) - qq--; - } - else - { - if (rr < yy / 2) - 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_centered_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_centered_quotient */ - return scm_i_bigint_centered_quotient (scm_i_long2big (xx), y); - } + return scm_integer_centered_quotient_iz (SCM_I_INUM (x), y); else if (SCM_REALP (y)) - return scm_i_inexact_centered_quotient (xx, SCM_REAL_VALUE (y)); + return scm_i_inexact_centered_quotient (SCM_I_INUM (x), + SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) return scm_i_exact_rational_centered_quotient (x, y); else @@ -2189,44 +2146,10 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-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_centered_quotient); - else if (SCM_UNLIKELY (yy == 1)) - return x; - else - { - SCM q = scm_i_mkbig (); - scm_t_inum rr; - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ - if (yy > 0) - { - rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), yy); - scm_remember_upto_here_1 (x); - if (rr < -yy / 2) - mpz_sub_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - } - else - { - rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), -yy); - scm_remember_upto_here_1 (x); - mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); - if (rr < yy / 2) - 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_centered_quotient_zi (x, SCM_I_INUM (y)); else if (SCM_BIGP (y)) - return scm_i_bigint_centered_quotient (x, y); + return scm_integer_centered_quotient_zz (x, y); else if (SCM_REALP (y)) return scm_i_inexact_centered_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); @@ -2276,49 +2199,6 @@ scm_i_inexact_centered_quotient (double x, double y) return scm_nan (); } -/* Assumes that both x and y are bigints, though - x might be able to fit into a fixnum. */ -static SCM -scm_i_bigint_centered_quotient (SCM x, SCM y) -{ - SCM q, r, min_r; - - /* 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 (); - - /* min_r will eventually become -abs(y)/2 */ - min_r = scm_i_mkbig (); - mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), - SCM_I_BIG_MPZ (y), 1); - - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ - if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) - { - mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - mpz_sub_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - } - else - { - mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - mpz_add_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - } - scm_remember_upto_here_2 (r, min_r); - return scm_i_normbig (q); -} - static SCM scm_i_exact_rational_centered_quotient (SCM x, SCM y) {