diff --git a/libguile/integers.c b/libguile/integers.c index 45f769f56..cca40bc8b 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -1742,3 +1742,78 @@ 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); } + +SCM +scm_integer_gcd_ii (scm_t_inum x, scm_t_inum y) +{ + scm_t_inum u = x < 0 ? -x : x; + scm_t_inum v = y < 0 ? -y : y; + scm_t_inum result; + if (x == 0) + result = v; + else if (y == 0) + result = u; + else + { + int k = 0; + /* Determine a common factor 2^k */ + while (((u | v) & 1) == 0) + { + k++; + u >>= 1; + v >>= 1; + } + /* Now, any factor 2^n can be eliminated */ + if ((u & 1) == 0) + while ((u & 1) == 0) + u >>= 1; + else + while ((v & 1) == 0) + v >>= 1; + /* Both u and v are now odd. Subtract the smaller one + from the larger one to produce an even number, remove + more factors of two, and repeat. */ + while (u != v) + { + if (u > v) + { + u -= v; + while ((u & 1) == 0) + u >>= 1; + } + else + { + v -= u; + while ((v & 1) == 0) + v >>= 1; + } + } + result = u << k; + } + return ulong_to_scm (result); +} + +SCM +scm_integer_gcd_zi (SCM x, scm_t_inum y) +{ + scm_t_bits result; + if (y == 0) + return scm_abs (x); + if (y < 0) + y = -y; + result = mpz_gcd_ui (NULL, SCM_I_BIG_MPZ (x), y); + scm_remember_upto_here_1 (x); + return ulong_to_scm (result); +} + +SCM +scm_integer_gcd_zz (SCM x, SCM y) +{ + mpz_t result, zx, zy; + mpz_init (result); + alias_bignum_to_mpz (scm_bignum (x), zx); + alias_bignum_to_mpz (scm_bignum (y), zy); + mpz_gcd (result, zx, zy); + scm_remember_upto_here_2 (x, y); + return take_mpz (result); +} diff --git a/libguile/integers.h b/libguile/integers.h index e1193c231..699525b3f 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -124,6 +124,10 @@ SCM_INTERNAL void scm_integer_round_divide_zi (SCM x, scm_t_inum y, SCM_INTERNAL void scm_integer_round_divide_zz (SCM x, SCM y, SCM *qp, SCM *rp); +SCM_INTERNAL SCM scm_integer_gcd_ii (scm_t_inum x, scm_t_inum y); +SCM_INTERNAL SCM scm_integer_gcd_zi (SCM x, scm_t_inum y); +SCM_INTERNAL SCM scm_integer_gcd_zz (SCM x, SCM y); + #endif /* SCM_INTEGERS_H */ diff --git a/libguile/numbers.c b/libguile/numbers.c index 32dad4d01..f0ddc9379 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -2782,68 +2782,15 @@ SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1, SCM scm_gcd (SCM x, SCM y) { - if (SCM_UNLIKELY (SCM_UNBNDP (y))) + if (SCM_UNBNDP (y)) return SCM_UNBNDP (x) ? SCM_INUM0 : scm_abs (x); - if (SCM_LIKELY (SCM_I_INUMP (x))) + if (SCM_I_INUMP (x)) { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum xx = SCM_I_INUM (x); - scm_t_inum yy = SCM_I_INUM (y); - scm_t_inum u = xx < 0 ? -xx : xx; - scm_t_inum v = yy < 0 ? -yy : yy; - scm_t_inum result; - if (SCM_UNLIKELY (xx == 0)) - result = v; - else if (SCM_UNLIKELY (yy == 0)) - result = u; - else - { - int k = 0; - /* Determine a common factor 2^k */ - while (((u | v) & 1) == 0) - { - k++; - u >>= 1; - v >>= 1; - } - /* Now, any factor 2^n can be eliminated */ - if ((u & 1) == 0) - while ((u & 1) == 0) - u >>= 1; - else - while ((v & 1) == 0) - v >>= 1; - /* Both u and v are now odd. Subtract the smaller one - from the larger one to produce an even number, remove - more factors of two, and repeat. */ - while (u != v) - { - if (u > v) - { - u -= v; - while ((u & 1) == 0) - u >>= 1; - } - else - { - v -= u; - while ((v & 1) == 0) - v >>= 1; - } - } - result = u << k; - } - return (SCM_POSFIXABLE (result) - ? SCM_I_MAKINUM (result) - : scm_i_inum2big (result)); - } + if (SCM_I_INUMP (y)) + return scm_integer_gcd_ii (SCM_I_INUM (x), SCM_I_INUM (y)); else if (SCM_BIGP (y)) - { - SCM_SWAP (x, y); - goto big_inum; - } + return scm_integer_gcd_zi (y, SCM_I_INUM (x)); else if (SCM_REALP (y) && scm_is_integer (y)) goto handle_inexacts; else @@ -2852,30 +2799,9 @@ scm_gcd (SCM x, SCM y) else if (SCM_BIGP (x)) { if (SCM_I_INUMP (y)) - { - scm_t_bits result; - scm_t_inum yy; - big_inum: - yy = SCM_I_INUM (y); - if (yy == 0) - return scm_abs (x); - if (yy < 0) - yy = -yy; - result = mpz_gcd_ui (NULL, SCM_I_BIG_MPZ (x), yy); - scm_remember_upto_here_1 (x); - return (SCM_POSFIXABLE (result) - ? SCM_I_MAKINUM (result) - : scm_from_unsigned_integer (result)); - } + return scm_integer_gcd_zi (x, SCM_I_INUM (y)); else if (SCM_BIGP (y)) - { - SCM result = scm_i_mkbig (); - mpz_gcd (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - return scm_i_normbig (result); - } + return scm_integer_gcd_zz (x, y); else if (SCM_REALP (y) && scm_is_integer (y)) goto handle_inexacts; else