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

Implement round-remainder with new integer lib

* libguile/integers.c (scm_integer_round_remainder_ii)
(scm_integer_round_remainder_iz, scm_integer_round_remainder_zi)
(scm_integer_round_remainder_zz): New internal functions.
(integer_round_remainder_zz): New helper.
(long_sign, bignum_cmp_long): New helpers.
* libguile/integers.h: Declare internal functions.
* libguile/numbers.c (scm_round_remainder): Use the new functions.
(scm_i_bigint_round_remainder): Remove unused helper.
This commit is contained in:
Andy Wingo 2021-12-13 11:53:34 +01:00
parent 9a358a9632
commit b114642640
3 changed files with 130 additions and 114 deletions

View file

@ -1503,3 +1503,119 @@ scm_integer_round_quotient_zz (SCM x, SCM y)
return scm_i_normbig (q);
}
static SCM
integer_round_remainder_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);
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_sub (r, r, zy);
scm_remember_upto_here_1 (y);
mpz_clear (q);
mpz_clear (r2);
return take_mpz (r);
}
SCM
scm_integer_round_remainder_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("round-remainder");
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)
r -= y;
else if (r2 <= -ay)
r += y;
}
else
{
if (r2 > ay)
r -= y;
else if (r2 < -ay)
r += y;
}
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_round_remainder_iz (scm_t_inum x, SCM y)
{
return integer_round_remainder_zz (long_to_bignum (x), scm_bignum (y));
}
SCM
scm_integer_round_remainder_zi (SCM x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("round-remainder");
mpz_t q, zx;
scm_t_inum r;
int needs_adjustment;
mpz_init (q);
alias_bignum_to_mpz (scm_bignum (x), zx);
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);
if (mpz_odd_p (q))
needs_adjustment = (2*r <= y);
else
needs_adjustment = (2*r < y);
}
scm_remember_upto_here_1 (x);
mpz_clear (q);
if (needs_adjustment)
r -= y;
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_round_remainder_zz (SCM x, SCM y)
{
return integer_round_remainder_zz (scm_bignum (x), scm_bignum (y));
}

View file

@ -110,6 +110,11 @@ 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);
SCM_INTERNAL SCM scm_integer_round_remainder_ii (scm_t_inum x, scm_t_inum y);
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);
#endif /* SCM_INTEGERS_H */

View file

@ -2533,7 +2533,6 @@ scm_i_exact_rational_round_quotient (SCM x, SCM y)
}
static SCM scm_i_inexact_round_remainder (double x, double y);
static SCM scm_i_bigint_round_remainder (SCM x, SCM y);
static SCM scm_i_exact_rational_round_remainder (SCM x, SCM y);
SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0,
@ -2556,53 +2555,15 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0,
"@end lisp")
#define FUNC_NAME s_scm_round_remainder
{
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_remainder);
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)
rr -= yy;
else if (r2 <= -ay)
rr += yy;
}
else
{
if (r2 > ay)
rr -= yy;
else if (r2 < -ay)
rr += yy;
}
return SCM_I_MAKINUM (rr);
}
}
if (SCM_I_INUMP (y))
return scm_integer_round_remainder_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_remainder */
return scm_i_bigint_round_remainder
(scm_i_long2big (xx), y);
}
return scm_integer_round_remainder_iz (SCM_I_INUM (x), y);
else if (SCM_REALP (y))
return scm_i_inexact_round_remainder (xx, SCM_REAL_VALUE (y));
return scm_i_inexact_round_remainder (SCM_I_INUM (x),
SCM_REAL_VALUE (y));
else if (SCM_FRACTIONP (y))
return scm_i_exact_rational_round_remainder (x, y);
else
@ -2611,43 +2572,10 @@ SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 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_remainder);
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);
if (mpz_odd_p (SCM_I_BIG_MPZ (q)))
needs_adjustment = (2*rr <= yy);
else
needs_adjustment = (2*rr < yy);
}
scm_remember_upto_here_2 (x, q);
if (needs_adjustment)
rr -= yy;
return SCM_I_MAKINUM (rr);
}
}
if (SCM_I_INUMP (y))
return scm_integer_round_remainder_zi (x, SCM_I_INUM (y));
else if (SCM_BIGP (y))
return scm_i_bigint_round_remainder (x, y);
return scm_integer_round_remainder_zz (x, y);
else if (SCM_REALP (y))
return scm_i_inexact_round_remainder
(scm_i_big2dbl (x), SCM_REAL_VALUE (y));
@ -2704,39 +2632,6 @@ scm_i_inexact_round_remainder (double x, double 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_remainder (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));
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);
scm_remember_upto_here_2 (q, r2);
if (needs_adjustment)
mpz_sub (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y));
scm_remember_upto_here_1 (y);
return scm_i_normbig (r);
}
static SCM
scm_i_exact_rational_round_remainder (SCM x, SCM y)
{