1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-05-20 11:40:18 +02:00

Implement centered-divide with new integer lib

* libguile/integers.c (scm_integer_centered_divide_ii)
(scm_integer_centered_divide_iz, scm_integer_centered_divide_zi)
(scm_integer_centered_divide_zz): New internal functions.
(integer_centered_divide_zz): New helper.
* libguile/integers.h: Declare internal functions.
* libguile/numbers.c (scm_centered_divide): Use the new functions.
(scm_i_bigint_centered_divide): Remove unused helper.
This commit is contained in:
Andy Wingo 2021-12-13 11:06:30 +01:00
parent 99b046d58b
commit ccb78fc7b1
3 changed files with 151 additions and 143 deletions

View file

@ -1231,3 +1231,136 @@ scm_integer_centered_remainder_zz (SCM x, SCM y)
{
return integer_centered_remainder_zz (scm_bignum (x), scm_bignum (y));
}
static void
integer_centered_divide_zz (struct scm_bignum *x, struct scm_bignum *y,
SCM *qp, SCM *rp)
{
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);
mpz_neg (min_r, min_r);
if (mpz_cmp (r, min_r) < 0)
{
mpz_sub_ui (q, q, 1);
mpz_add (r, r, zy);
}
}
else
{
mpz_fdiv_qr (q, r, zx, zy);
if (mpz_cmp (r, min_r) < 0)
{
mpz_add_ui (q, q, 1);
mpz_sub (r, r, zy);
}
}
scm_remember_upto_here_2 (x, y);
mpz_clear (min_r);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
void
scm_integer_centered_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("centered-divide");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
if (x > 0)
{
if (y > 0)
{
if (r >= (y + 1) / 2)
{ q++; r -= y; }
}
else
{
if (r >= (1 - y) / 2)
{ q--; r += y; }
}
}
else
{
if (y > 0)
{
if (r < -y / 2)
{ q--; r += y; }
}
else
{
if (r < y / 2)
{ q++; r -= y; }
}
}
*qp = long_to_scm (q);
*rp = SCM_I_MAKINUM (r);
}
void
scm_integer_centered_divide_iz (scm_t_inum x, SCM y, SCM *qp, SCM *rp)
{
integer_centered_divide_zz (long_to_bignum (x), scm_bignum (y), qp, rp);
}
void
scm_integer_centered_divide_zi (SCM x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("centered-divide");
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);
if (r < -y / 2)
{
mpz_sub_ui (q, q, 1);
r += y;
}
}
else
{
r = - mpz_cdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
if (r < y / 2)
{
mpz_add_ui (q, q, 1);
r -= y;
}
}
scm_remember_upto_here_1 (x);
*qp = take_mpz (q);
*rp = SCM_I_MAKINUM (r);
}
void
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);
}

View file

@ -96,6 +96,15 @@ SCM_INTERNAL SCM scm_integer_centered_remainder_iz (scm_t_inum x, SCM y);
SCM_INTERNAL SCM scm_integer_centered_remainder_zi (SCM x, scm_t_inum y);
SCM_INTERNAL SCM scm_integer_centered_remainder_zz (SCM x, SCM y);
SCM_INTERNAL void scm_integer_centered_divide_ii (scm_t_inum x, scm_t_inum y,
SCM *qp, SCM *rp);
SCM_INTERNAL void scm_integer_centered_divide_iz (scm_t_inum x, SCM y,
SCM *qp, SCM *rp);
SCM_INTERNAL void scm_integer_centered_divide_zi (SCM x, scm_t_inum y,
SCM *qp, SCM *rp);
SCM_INTERNAL void scm_integer_centered_divide_zz (SCM x, SCM y,
SCM *qp, SCM *rp);
#endif /* SCM_INTEGERS_H */

View file

@ -2320,7 +2320,6 @@ scm_i_exact_rational_centered_remainder (SCM x, SCM y)
static void scm_i_inexact_centered_divide (double x, double y,
SCM *qp, SCM *rp);
static void scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp);
static void scm_i_exact_rational_centered_divide (SCM x, SCM y,
SCM *qp, SCM *rp);
@ -2352,57 +2351,15 @@ SCM_PRIMITIVE_GENERIC (scm_i_centered_divide, "centered/", 2, 0, 0,
void
scm_centered_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_centered_divide);
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++; rr -= yy; }
}
else
{
if (rr >= (1 - yy) / 2)
{ qq--; rr += yy; }
}
}
else
{
if (SCM_LIKELY (yy > 0))
{
if (rr < -yy / 2)
{ qq--; rr += yy; }
}
else
{
if (rr < yy / 2)
{ 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_centered_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_centered_divide */
scm_i_bigint_centered_divide (scm_i_long2big (xx), y, qp, rp);
scm_integer_centered_divide_iz (SCM_I_INUM (x), y, qp, rp);
else if (SCM_REALP (y))
scm_i_inexact_centered_divide (xx, SCM_REAL_VALUE (y), qp, rp);
scm_i_inexact_centered_divide (SCM_I_INUM (x), SCM_REAL_VALUE (y),
qp, rp);
else if (SCM_FRACTIONP (y))
scm_i_exact_rational_centered_divide (x, y, qp, rp);
else
@ -2411,49 +2368,10 @@ scm_centered_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_centered_divide);
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);
rr += yy;
}
}
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);
rr -= yy;
}
}
*qp = scm_i_normbig (q);
*rp = SCM_I_MAKINUM (rr);
}
}
if (SCM_I_INUMP (y))
scm_integer_centered_divide_zi (x, SCM_I_INUM (y), qp, rp);
else if (SCM_BIGP (y))
scm_i_bigint_centered_divide (x, y, qp, rp);
scm_integer_centered_divide_zz (x, y, qp, rp);
else if (SCM_REALP (y))
scm_i_inexact_centered_divide (scm_i_big2dbl (x), SCM_REAL_VALUE (y),
qp, rp);
@ -2507,58 +2425,6 @@ scm_i_inexact_centered_divide (double x, double y, SCM *qp, SCM *rp)
*rp = scm_i_from_double (r);
}
/* Assumes that both x and y are bigints, though
x might be able to fit into a fixnum. */
static void
scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{
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));
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);
mpz_add (SCM_I_BIG_MPZ (r),
SCM_I_BIG_MPZ (r),
SCM_I_BIG_MPZ (y));
}
}
else
{
mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (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);
mpz_sub (SCM_I_BIG_MPZ (r),
SCM_I_BIG_MPZ (r),
SCM_I_BIG_MPZ (y));
}
}
scm_remember_upto_here_2 (x, y);
*qp = scm_i_normbig (q);
*rp = scm_i_normbig (r);
}
static void
scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp)
{