summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAndy Wingo <wingo@pobox.com>2021-12-13 12:12:54 +0100
committerAndy Wingo <wingo@pobox.com>2022-01-13 09:37:16 +0100
commit025c7c80451c9394cff1e15ce278d34a01a98022 (patch)
tree49f28ee0d5ddf1a184ff3a30df365d8f36f26fc8
parentb114642640a6b4a34de243cb6a4681ac3e140bb4 (diff)
downloadguile-025c7c80451c9394cff1e15ce278d34a01a98022.tar.gz
Implement round-divide with new integer lib
* libguile/integers.c (scm_integer_round_divide_ii) (scm_integer_round_divide_iz, scm_integer_round_divide_zi) (scm_integer_round_divide_zz): New internal functions. (integer_round_divide_zz): New helper. * libguile/integers.h: Declare internal functions. * libguile/numbers.c (scm_round_divide): Use the new functions. (scm_i_bigint_round_divide): Remove unused helper.
-rw-r--r--libguile/integers.c123
-rw-r--r--libguile/integers.h9
-rw-r--r--libguile/numbers.c131
3 files changed, 140 insertions, 123 deletions
diff --git a/libguile/integers.c b/libguile/integers.c
index 715d2ed6f..45f769f56 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -1619,3 +1619,126 @@ scm_integer_round_remainder_zz (SCM x, SCM y)
{
return integer_round_remainder_zz (scm_bignum (x), scm_bignum (y));
}
+
+static void
+integer_round_divide_zz (struct scm_bignum *x, struct scm_bignum *y,
+ SCM *qp, SCM *rp)
+{
+ 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);
+ scm_remember_upto_here_1 (x);
+ mpz_mul_2exp (r2, r, 1); /* r2 = 2*r */
+
+ cmp = mpz_cmpabs (r2, zy);
+ if (mpz_odd_p (q))
+ needs_adjustment = (cmp >= 0);
+ else
+ needs_adjustment = (cmp > 0);
+
+ if (needs_adjustment)
+ {
+ mpz_add_ui (q, q, 1);
+ mpz_sub (r, r, zy);
+ }
+
+ scm_remember_upto_here_1 (y);
+ mpz_clear (r2);
+ *qp = take_mpz (q);
+ *rp = take_mpz (r);
+}
+
+void
+scm_integer_round_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
+{
+ if (y == 0)
+ scm_num_overflow ("round-divide");
+
+ 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++; r -= y; }
+ else if (r2 <= -ay)
+ { q--; r += y; }
+ }
+ else
+ {
+ if (r2 > ay)
+ { q++; r -= y; }
+ else if (r2 < -ay)
+ { q--; r += y; }
+ }
+ *qp = long_to_scm (q);
+ *rp = SCM_I_MAKINUM (r);
+}
+
+void
+scm_integer_round_divide_iz (scm_t_inum x, SCM y, SCM *qp, SCM *rp)
+{
+ integer_round_divide_zz (long_to_bignum (x), scm_bignum (y), qp, rp);
+}
+
+void
+scm_integer_round_divide_zi (SCM x, scm_t_inum y, SCM *qp, SCM *rp)
+{
+ if (y == 0)
+ scm_num_overflow ("round-divide");
+
+ 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);
+ r -= y;
+ }
+ *qp = take_mpz (q);
+ *rp = SCM_I_MAKINUM (r);
+}
+
+void
+scm_integer_round_divide_zz (SCM x, SCM y, SCM *qp, SCM *rp)
+{
+ integer_round_divide_zz (scm_bignum (x), scm_bignum (y), qp, rp);
+}
diff --git a/libguile/integers.h b/libguile/integers.h
index 6628b444d..e1193c231 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -115,6 +115,15 @@ SCM_INTERNAL SCM scm_integer_round_remainder_iz (scm_t_inum x, SCM y);
SCM_INTERNAL SCM scm_integer_round_remainder_zi (SCM x, scm_t_inum y);
SCM_INTERNAL SCM scm_integer_round_remainder_zz (SCM x, SCM y);
+SCM_INTERNAL void scm_integer_round_divide_ii (scm_t_inum x, scm_t_inum y,
+ SCM *qp, SCM *rp);
+SCM_INTERNAL void scm_integer_round_divide_iz (scm_t_inum x, SCM y,
+ SCM *qp, SCM *rp);
+SCM_INTERNAL void scm_integer_round_divide_zi (SCM x, scm_t_inum y,
+ SCM *qp, SCM *rp);
+SCM_INTERNAL void scm_integer_round_divide_zz (SCM x, SCM y,
+ SCM *qp, SCM *rp);
+
#endif /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 55dfd58cf..32dad4d01 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -2644,7 +2644,6 @@ scm_i_exact_rational_round_remainder (SCM x, SCM y)
static void scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp);
-static void scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp);
static void scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp);
SCM_PRIMITIVE_GENERIC (scm_i_round_divide, "round/", 2, 0, 0,
@@ -2679,54 +2678,14 @@ SCM_PRIMITIVE_GENERIC (scm_i_round_divide, "round/", 2, 0, 0,
void
scm_round_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{
- 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_divide);
- 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++; rr -= yy; }
- else if (r2 <= -ay)
- { qq--; rr += yy; }
- }
- else
- {
- if (r2 > ay)
- { qq++; rr -= yy; }
- else if (r2 < -ay)
- { qq--; rr += yy; }
- }
- if (SCM_LIKELY (SCM_FIXABLE (qq)))
- *qp = SCM_I_MAKINUM (qq);
- else
- *qp = scm_i_inum2big (qq);
- *rp = SCM_I_MAKINUM (rr);
- }
- }
+ if (SCM_I_INUMP (y))
+ scm_integer_round_divide_ii (SCM_I_INUM (x), SCM_I_INUM (y), qp, rp);
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_divide */
- scm_i_bigint_round_divide (scm_i_long2big (SCM_I_INUM (x)), y, qp, rp);
+ scm_integer_round_divide_iz (SCM_I_INUM (x), y, qp, rp);
else if (SCM_REALP (y))
- scm_i_inexact_round_divide (xx, SCM_REAL_VALUE (y), qp, rp);
+ scm_i_inexact_round_divide (SCM_I_INUM (x), SCM_REAL_VALUE (y), qp, rp);
else if (SCM_FRACTIONP (y))
scm_i_exact_rational_round_divide (x, y, qp, rp);
else
@@ -2735,48 +2694,10 @@ scm_round_divide (SCM x, SCM y, SCM *qp, SCM *rp)
}
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_divide);
- 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);
- rr -= yy;
- }
- *qp = scm_i_normbig (q);
- *rp = SCM_I_MAKINUM (rr);
- }
- }
+ if (SCM_I_INUMP (y))
+ scm_integer_round_divide_zi (x, SCM_I_INUM (y), qp, rp);
else if (SCM_BIGP (y))
- scm_i_bigint_round_divide (x, y, qp, rp);
+ scm_integer_round_divide_zz (x, y, qp, rp);
else if (SCM_REALP (y))
scm_i_inexact_round_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y),
qp, rp);
@@ -2826,42 +2747,6 @@ scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp)
}
}
-/* Assumes that both x and y are bigints, though
- x might be able to fit into a fixnum. */
-static void
-scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp)
-{
- 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));
- scm_remember_upto_here_1 (x);
- mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*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);
-
- if (needs_adjustment)
- {
- mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
- mpz_sub (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y));
- }
-
- scm_remember_upto_here_2 (r2, y);
- *qp = scm_i_normbig (q);
- *rp = scm_i_normbig (r);
-}
-
static void
scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{