diff options
author | Andy Wingo <wingo@pobox.com> | 2022-01-04 22:32:27 +0100 |
---|---|---|
committer | Andy Wingo <wingo@pobox.com> | 2022-01-13 09:37:16 +0100 |
commit | 8b2d58b993f1aa7625156d3a5e66efd9acef7262 (patch) | |
tree | 249dc2c538952f8b223b7087edc7d1be054375d9 | |
parent | 3e08c9cec03660048cbc4473b1e9ed421cd01918 (diff) | |
download | guile-8b2d58b993f1aa7625156d3a5e66efd9acef7262.tar.gz |
Clean up scm_divide
* libguile/integers.h:
* libguile/integers.c (scm_is_integer_divisible_ii):
(scm_is_integer_divisible_zi):
(scm_is_integer_divisible_zz): New helpers.
* libguile/numbers.c (invert, divide, complex_div): New helpers for
scm_divide.
(scm_divide): Adapt.
-rw-r--r-- | libguile/integers.c | 85 | ||||
-rw-r--r-- | libguile/integers.h | 12 | ||||
-rw-r--r-- | libguile/numbers.c | 378 |
3 files changed, 252 insertions, 223 deletions
diff --git a/libguile/integers.c b/libguile/integers.c index 715c210f5..4b26ea3e0 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -2658,3 +2658,88 @@ scm_integer_mul_zz (struct scm_bignum *x, struct scm_bignum *y) scm_remember_upto_here_2 (x, y); return take_mpz (result); } + +int +scm_is_integer_divisible_ii (scm_t_inum x, scm_t_inum y) +{ + ASSERT (y != 0); + return (x % y) == 0; +} + +int +scm_is_integer_divisible_zi (struct scm_bignum *x, scm_t_inum y) +{ + ASSERT (y != 0); + switch (y) + { + case -1: + case 1: + return 1; + default: + { + scm_t_inum abs_y = y < 0 ? -y : y; + mpz_t zx; + alias_bignum_to_mpz (x, zx); + int divisible = mpz_divisible_ui_p (zx, abs_y); + scm_remember_upto_here_1 (x); + return divisible; + } + } +} + +int +scm_is_integer_divisible_zz (struct scm_bignum *x, struct scm_bignum *y) +{ + mpz_t zx, zy; + alias_bignum_to_mpz (x, zx); + alias_bignum_to_mpz (y, zy); + int divisible_p = mpz_divisible_p (zx, zy); + scm_remember_upto_here_2 (x, y); + return divisible_p; +} + +SCM +scm_integer_exact_quotient_ii (scm_t_inum n, scm_t_inum d) +{ + return scm_integer_truncate_quotient_ii (n, d); +} + +/* Return the exact integer q such that n = q*d, for exact integers n + and d, where d is known in advance to divide n evenly (with zero + remainder). For large integers, this can be computed more + efficiently than when the remainder is unknown. */ +SCM +scm_integer_exact_quotient_zi (struct scm_bignum *n, scm_t_inum d) +{ + if (SCM_UNLIKELY (d == 0)) + scm_num_overflow ("quotient"); + else if (SCM_UNLIKELY (d == 1)) + return scm_from_bignum (n); + + mpz_t q, zn; + mpz_init (q); + alias_bignum_to_mpz (n, zn); + if (d > 0) + mpz_divexact_ui (q, zn, d); + else + { + mpz_divexact_ui (q, zn, -d); + mpz_neg (q, q); + } + scm_remember_upto_here_1 (n); + return take_mpz (q); +} + +SCM +scm_integer_exact_quotient_zz (struct scm_bignum *n, struct scm_bignum *d) +{ + mpz_t q, zn, zd; + mpz_init (q); + alias_bignum_to_mpz (n, zn); + alias_bignum_to_mpz (d, zd); + + mpz_divexact (q, zn, zd); + scm_remember_upto_here_2 (n, d); + return take_mpz (q); +} + diff --git a/libguile/integers.h b/libguile/integers.h index 8795a6288..98c19b90f 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -184,6 +184,18 @@ SCM_INTERNAL SCM scm_integer_mul_ii (scm_t_inum x, scm_t_inum y); SCM_INTERNAL SCM scm_integer_mul_zi (struct scm_bignum *x, scm_t_inum y); SCM_INTERNAL SCM scm_integer_mul_zz (struct scm_bignum *x, struct scm_bignum *y); +SCM_INTERNAL int scm_is_integer_divisible_ii (scm_t_inum x, scm_t_inum y); +SCM_INTERNAL int scm_is_integer_divisible_zi (struct scm_bignum *x, + scm_t_inum y); +SCM_INTERNAL int scm_is_integer_divisible_zz (struct scm_bignum *x, + struct scm_bignum *y); + +SCM_INTERNAL SCM scm_integer_exact_quotient_ii (scm_t_inum n, scm_t_inum d); +SCM_INTERNAL SCM scm_integer_exact_quotient_zi (struct scm_bignum *n, + scm_t_inum d); +SCM_INTERNAL SCM scm_integer_exact_quotient_zz (struct scm_bignum *n, + struct scm_bignum *d); + #endif /* SCM_INTEGERS_H */ diff --git a/libguile/numbers.c b/libguile/numbers.c index c544b08d2..17131b9ba 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -5576,259 +5576,163 @@ arising out of or in connection with the use or performance of this software. ****************************************************************/ -SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1, - (SCM x, SCM y, SCM rest), - "Divide the first argument by the product of the remaining\n" - "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n" - "returned.") -#define FUNC_NAME s_scm_i_divide +static SCM +invert (SCM x) { - while (!scm_is_null (rest)) - { x = scm_divide (x, y); - y = scm_car (rest); - rest = scm_cdr (rest); + if (SCM_I_INUMP (x)) + switch (SCM_I_INUM (x)) + { + case -1: return x; + case 0: scm_num_overflow ("divide"); + case 1: return x; + default: return scm_i_make_ratio_already_reduced (SCM_INUM1, x); + } + else if (SCM_BIGP (x)) + return scm_i_make_ratio_already_reduced (SCM_INUM1, x); + else if (SCM_REALP (x)) + return scm_i_from_double (1.0 / SCM_REAL_VALUE (x)); + else if (SCM_COMPLEXP (x)) + { + double r = SCM_COMPLEX_REAL (x); + double i = SCM_COMPLEX_IMAG (x); + if (fabs(r) <= fabs(i)) + { + double t = r / i; + double d = i * (1.0 + t * t); + return scm_c_make_rectangular (t / d, -1.0 / d); + } + else + { + double t = i / r; + double d = r * (1.0 + t * t); + return scm_c_make_rectangular (1.0 / d, -t / d); + } } - return scm_divide (x, y); + else if (SCM_FRACTIONP (x)) + return scm_i_make_ratio_already_reduced (SCM_FRACTION_DENOMINATOR (x), + SCM_FRACTION_NUMERATOR (x)); + else + abort (); /* Unreachable. */ } -#undef FUNC_NAME - -#define s_divide s_scm_i_divide -#define g_divide g_scm_i_divide -SCM -scm_divide (SCM x, SCM y) -#define FUNC_NAME s_divide +static SCM +complex_div (double a, SCM y) { - double a; - - if (SCM_UNLIKELY (SCM_UNBNDP (y))) + double r = SCM_COMPLEX_REAL (y); + double i = SCM_COMPLEX_IMAG (y); + if (fabs(r) <= fabs(i)) { - if (SCM_UNBNDP (x)) - return scm_wta_dispatch_0 (g_divide, s_divide); - else if (SCM_I_INUMP (x)) - { - scm_t_inum xx = SCM_I_INUM (x); - if (xx == 1 || xx == -1) - return x; - else if (xx == 0) - scm_num_overflow (s_divide); - else - return scm_i_make_ratio_already_reduced (SCM_INUM1, x); - } - else if (SCM_BIGP (x)) - return scm_i_make_ratio_already_reduced (SCM_INUM1, x); - else if (SCM_REALP (x)) - return scm_i_from_double (1.0 / SCM_REAL_VALUE (x)); - else if (SCM_COMPLEXP (x)) - { - double r = SCM_COMPLEX_REAL (x); - double i = SCM_COMPLEX_IMAG (x); - if (fabs(r) <= fabs(i)) - { - double t = r / i; - double d = i * (1.0 + t * t); - return scm_c_make_rectangular (t / d, -1.0 / d); - } - else - { - double t = i / r; - double d = r * (1.0 + t * t); - return scm_c_make_rectangular (1.0 / d, -t / d); - } - } - else if (SCM_FRACTIONP (x)) - return scm_i_make_ratio_already_reduced (SCM_FRACTION_DENOMINATOR (x), - SCM_FRACTION_NUMERATOR (x)); - else - return scm_wta_dispatch_1 (g_divide, x, SCM_ARG1, s_divide); + double t = r / i; + double d = i * (1.0 + t * t); + return scm_c_make_rectangular ((a * t) / d, -a / d); } + else + { + double t = i / r; + double d = r * (1.0 + t * t); + return scm_c_make_rectangular (a / d, -(a * t) / d); + } +} - if (SCM_LIKELY (SCM_I_INUMP (x))) +static SCM +divide (SCM x, SCM y) +{ + if (scm_is_eq (y, SCM_INUM0)) + scm_num_overflow ("divide"); + + 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 (yy == 0) - scm_num_overflow (s_divide); - else if (xx % yy != 0) - return scm_i_make_ratio (x, y); - else - { - scm_t_inum z = xx / yy; - if (SCM_FIXABLE (z)) - return SCM_I_MAKINUM (z); - else - return scm_i_inum2big (z); - } - } + if (scm_is_eq (x, SCM_INUM1)) + return invert (y); + if (SCM_I_INUMP (y)) + return scm_is_integer_divisible_ii (SCM_I_INUM (x), SCM_I_INUM (y)) + ? scm_integer_exact_quotient_ii (SCM_I_INUM (x), SCM_I_INUM (y)) + : scm_i_make_ratio (x, y); else if (SCM_BIGP (y)) return scm_i_make_ratio (x, y); else if (SCM_REALP (y)) /* FIXME: Precision may be lost here due to: (1) The cast from 'scm_t_inum' to 'double' (2) Double rounding */ - return scm_i_from_double ((double) xx / SCM_REAL_VALUE (y)); + return scm_i_from_double ((double) SCM_I_INUM (x) / SCM_REAL_VALUE (y)); else if (SCM_COMPLEXP (y)) - { - a = xx; - complex_div: /* y _must_ be a complex number */ - { - double r = SCM_COMPLEX_REAL (y); - double i = SCM_COMPLEX_IMAG (y); - if (fabs(r) <= fabs(i)) - { - double t = r / i; - double d = i * (1.0 + t * t); - return scm_c_make_rectangular ((a * t) / d, -a / d); - } - else - { - double t = i / r; - double d = r * (1.0 + t * t); - return scm_c_make_rectangular (a / d, -(a * t) / d); - } - } - } + return complex_div (SCM_I_INUM (x), y); else if (SCM_FRACTIONP (y)) /* a / b/c = ac / b */ return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), SCM_FRACTION_NUMERATOR (y)); else - return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); + abort (); /* Unreachable. */ } else if (SCM_BIGP (x)) { if (SCM_I_INUMP (y)) - { - scm_t_inum yy = SCM_I_INUM (y); - if (yy == 0) - scm_num_overflow (s_divide); - else if (yy == 1) - return x; - else - { - /* FIXME: HMM, what are the relative performance issues here? - We need to test. Is it faster on average to test - divisible_p, then perform whichever operation, or is it - faster to perform the integer div opportunistically and - switch to real if there's a remainder? For now we take the - middle ground: test, then if divisible, use the faster div - func. */ - - scm_t_inum abs_yy = yy < 0 ? -yy : yy; - int divisible_p = mpz_divisible_ui_p (SCM_I_BIG_MPZ (x), abs_yy); - - if (divisible_p) - { - SCM result = scm_i_mkbig (); - mpz_divexact_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), abs_yy); - scm_remember_upto_here_1 (x); - if (yy < 0) - mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result)); - return scm_i_normbig (result); - } - else - return scm_i_make_ratio (x, y); - } - } + return scm_is_integer_divisible_zi (scm_bignum (x), SCM_I_INUM (y)) + ? scm_integer_exact_quotient_zi (scm_bignum (x), SCM_I_INUM (y)) + : scm_i_make_ratio (x, y); else if (SCM_BIGP (y)) - { - int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - if (divisible_p) - { - SCM result = scm_i_mkbig (); - mpz_divexact (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); - } - else - return scm_i_make_ratio (x, y); - } + return scm_is_integer_divisible_zz (scm_bignum (x), scm_bignum (y)) + ? scm_integer_exact_quotient_zz (scm_bignum (x), scm_bignum (y)) + : scm_i_make_ratio (x, y); else if (SCM_REALP (y)) /* FIXME: Precision may be lost here due to: - (1) scm_i_big2dbl (2) Double rounding */ - return scm_i_from_double (scm_i_big2dbl (x) / SCM_REAL_VALUE (y)); + (1) scm_integer_to_double_z (2) Double rounding */ + return scm_i_from_double (scm_integer_to_double_z (scm_bignum (x)) + / SCM_REAL_VALUE (y)); else if (SCM_COMPLEXP (y)) - { - a = scm_i_big2dbl (x); - goto complex_div; - } + return complex_div (scm_integer_to_double_z (scm_bignum (x)), y); else if (SCM_FRACTIONP (y)) return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), SCM_FRACTION_NUMERATOR (y)); else - return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); + abort (); /* Unreachable. */ } else if (SCM_REALP (x)) { double rx = SCM_REAL_VALUE (x); if (SCM_I_INUMP (y)) - { - scm_t_inum yy = SCM_I_INUM (y); - if (yy == 0) - scm_num_overflow (s_divide); - else - /* FIXME: Precision may be lost here due to: - (1) The cast from 'scm_t_inum' to 'double' - (2) Double rounding */ - return scm_i_from_double (rx / (double) yy); - } + /* FIXME: Precision may be lost here due to: + (1) The cast from 'scm_t_inum' to 'double' + (2) Double rounding */ + return scm_i_from_double (rx / (double) SCM_I_INUM (y)); else if (SCM_BIGP (y)) - { - /* FIXME: Precision may be lost here due to: - (1) The conversion from bignum to double - (2) Double rounding */ - double dby = mpz_get_d (SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_1 (y); - return scm_i_from_double (rx / dby); - } + /* FIXME: Precision may be lost here due to: + (1) The conversion from bignum to double + (2) Double rounding */ + return scm_i_from_double (rx / scm_integer_to_double_z (scm_bignum (y))); else if (SCM_REALP (y)) return scm_i_from_double (rx / SCM_REAL_VALUE (y)); else if (SCM_COMPLEXP (y)) - { - a = rx; - goto complex_div; - } + return complex_div (rx, y); else if (SCM_FRACTIONP (y)) return scm_i_from_double (rx / scm_i_fraction2double (y)); else - return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); + abort () ; /* Unreachable. */ } else if (SCM_COMPLEXP (x)) { double rx = SCM_COMPLEX_REAL (x); double ix = SCM_COMPLEX_IMAG (x); if (SCM_I_INUMP (y)) - { - scm_t_inum yy = SCM_I_INUM (y); - if (yy == 0) - scm_num_overflow (s_divide); - else - { - /* FIXME: Precision may be lost here due to: - (1) The conversion from 'scm_t_inum' to double - (2) Double rounding */ - double d = yy; - return scm_c_make_rectangular (rx / d, ix / d); - } + { + /* FIXME: Precision may be lost here due to: + (1) The conversion from 'scm_t_inum' to double + (2) Double rounding */ + double d = SCM_I_INUM (y); + return scm_c_make_rectangular (rx / d, ix / d); } else if (SCM_BIGP (y)) { /* FIXME: Precision may be lost here due to: (1) The conversion from bignum to double (2) Double rounding */ - double dby = mpz_get_d (SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_1 (y); - return scm_c_make_rectangular (rx / dby, ix / dby); + double d = scm_integer_to_double_z (scm_bignum (y)); + return scm_c_make_rectangular (rx / d, ix / d); } else if (SCM_REALP (y)) { - double yy = SCM_REAL_VALUE (y); - return scm_c_make_rectangular (rx / yy, ix / yy); + double d = SCM_REAL_VALUE (y); + return scm_c_make_rectangular (rx / d, ix / d); } else if (SCM_COMPLEXP (y)) { @@ -5838,13 +5742,15 @@ scm_divide (SCM x, SCM y) { double t = ry / iy; double d = iy * (1.0 + t * t); - return scm_c_make_rectangular ((rx * t + ix) / d, (ix * t - rx) / d); + return scm_c_make_rectangular ((rx * t + ix) / d, + (ix * t - rx) / d); } else { double t = iy / ry; double d = ry * (1.0 + t * t); - return scm_c_make_rectangular ((rx + ix * t) / d, (ix - rx * t) / d); + return scm_c_make_rectangular ((rx + ix * t) / d, + (ix - rx * t) / d); } } else if (SCM_FRACTIONP (y)) @@ -5852,28 +5758,17 @@ scm_divide (SCM x, SCM y) /* FIXME: Precision may be lost here due to: (1) The conversion from fraction to double (2) Double rounding */ - double yy = scm_i_fraction2double (y); - return scm_c_make_rectangular (rx / yy, ix / yy); + double d = scm_i_fraction2double (y); + return scm_c_make_rectangular (rx / d, ix / d); } else - return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); + abort (); /* Unreachable. */ } else if (SCM_FRACTIONP (x)) { - if (SCM_I_INUMP (y)) - { - scm_t_inum yy = SCM_I_INUM (y); - if (yy == 0) - scm_num_overflow (s_divide); - else - return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x), - scm_product (SCM_FRACTION_DENOMINATOR (x), y)); - } - else if (SCM_BIGP (y)) - { - return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x), - scm_product (SCM_FRACTION_DENOMINATOR (x), y)); - } + if (scm_is_exact_integer (y)) + return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x), + scm_product (SCM_FRACTION_DENOMINATOR (x), y)); else if (SCM_REALP (y)) /* FIXME: Precision may be lost here due to: (1) The conversion from fraction to double @@ -5881,24 +5776,61 @@ scm_divide (SCM x, SCM y) return scm_i_from_double (scm_i_fraction2double (x) / SCM_REAL_VALUE (y)); else if (SCM_COMPLEXP (y)) - { - /* FIXME: Precision may be lost here due to: - (1) The conversion from fraction to double - (2) Double rounding */ - a = scm_i_fraction2double (x); - goto complex_div; - } + /* FIXME: Precision may be lost here due to: + (1) The conversion from fraction to double + (2) Double rounding */ + return complex_div (scm_i_fraction2double (x), y); else if (SCM_FRACTIONP (y)) - return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)), - scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))); + return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (y)), + scm_product (SCM_FRACTION_NUMERATOR (y), + SCM_FRACTION_DENOMINATOR (x))); else - return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); + abort (); /* Unreachable. */ } else - return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARG1, s_divide); + abort (); /* Unreachable. */ +} + +SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1, + (SCM x, SCM y, SCM rest), + "Divide the first argument by the product of the remaining\n" + "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n" + "returned.") +#define FUNC_NAME s_scm_i_divide +{ + while (!scm_is_null (rest)) + { x = scm_divide (x, y); + y = scm_car (rest); + rest = scm_cdr (rest); + } + return scm_divide (x, y); } #undef FUNC_NAME + +SCM +scm_divide (SCM x, SCM y) +{ + if (SCM_UNBNDP (y)) + { + if (SCM_UNBNDP (x)) + return scm_wta_dispatch_0 (g_scm_i_divide, s_scm_i_divide); + if (SCM_NUMBERP (x)) + return invert (x); + else + return scm_wta_dispatch_1 (g_scm_i_divide, x, SCM_ARG1, + s_scm_i_divide); + } + if (!SCM_NUMBERP (x)) + return scm_wta_dispatch_2 (g_scm_i_divide, x, y, SCM_ARG1, + s_scm_i_divide); + if (!SCM_NUMBERP (y)) + return scm_wta_dispatch_2 (g_scm_i_divide, x, y, SCM_ARG2, + s_scm_i_divide); + + return divide (x, y); +} double scm_c_truncate (double x) |