summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAndy Wingo <wingo@pobox.com>2021-12-13 11:42:17 +0100
committerAndy Wingo <wingo@pobox.com>2022-01-13 09:37:16 +0100
commit9a358a9632f6ef9c58f4d080389cd8bd4c9b0edc (patch)
treeb8389e163dd1020849e221721084b30c829567f4
parentccb78fc7b171d1f2bf1f6eeb9e4cffcea318a838 (diff)
downloadguile-9a358a9632f6ef9c58f4d080389cd8bd4c9b0edc.tar.gz
Implement round-quotient with new integer lib
* libguile/integers.c (scm_integer_round_quotient_ii) (scm_integer_round_quotient_iz, scm_integer_round_quotient_zi) (scm_integer_round_quotient_zz): New internal functions. (integer_round_quotient_zz): New helper. (long_sign, bignum_cmp_long): New helpers. * libguile/integers.h: Declare internal functions. * libguile/numbers.c (scm_round_quotient): Use the new functions. (scm_i_bigint_round_quotient): Remove unused helper.
-rw-r--r--libguile/integers.c141
-rw-r--r--libguile/integers.h5
-rw-r--r--libguile/numbers.c127
3 files changed, 154 insertions, 119 deletions
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)
{