From 807de5ec6dab40fc699f7e3353de36ae1ee36604 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Tue, 1 Mar 2011 16:03:59 -0500 Subject: [PATCH 2/7] Optimize and simplify fractions code * libguile/numbers.c (scm_exact_integer_quotient): New internal static function that computes the quotient of two exact integers when the remainder is known in advance to be zero. For large integers this can be implemented more efficiently than when the remainder is unknown. (scm_i_make_ratio_already_reduced): New internal static function that creates a ratio when the numerator and denominator are already known to be reduced to lowest terms (i.e. when their gcd is 1). This can be used in several places to avoid unnecessary gcd computations. (scm_i_make_ratio): Rewrite in terms of scm_i_make_ratio_already_reduced. Don't bother checking to see if the denominator divides the numerator evenly. This is wasted effort in the common case. Instead, compute the gcd, reduce to lowest terms (using scm_exact_integer_quotient), and let scm_i_make_ratio_already_reduced do the integer check (by checking for a unit denominator). (scm_integer_expt): Optimize fraction case by (a/b)^n ==> (a^n)/(b^n), to avoid needless effort to reduce intermediate products to lowest terms. If a and b have no common factors, it follows that a^n and b^n don't have common factors. Use scm_i_make_ratio_already_reduced to construct the final result, so that no gcd computations are needed to exponentiate a fraction. (scm_abs, scm_magnitude, scm_difference): Use scm_i_make_ratio_already_reduced to avoid wasteful gcd computations when negating fractions. (do_divide): Use scm_i_make_ratio_already_reduced to avoid wasteful gcd computations when taking the reciprocal of a fraction or integer. When dividing exact integers, don't bother checking to see if the denominator divides the numerator evenly. Instead, let scm_i_make_ratio do the job. * test-suite/tests/numbers.test (expt, integer-expt): Add tests. --- libguile/numbers.c | 287 ++++++++++++++++++++++------------------- test-suite/tests/numbers.test | 6 + 2 files changed, 161 insertions(+), 132 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index aa768b4..80b46df 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -413,96 +413,58 @@ scm_i_mpz2num (mpz_t b) /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */ static SCM scm_divide2real (SCM x, SCM y); +/* scm_i_make_ratio_already_reduced makes a ratio more efficiently, + when the following conditions are known in advance: + + 1. numerator and denominator are exact integers + 2. numerator and denominator are reduced to lowest terms: gcd(n,d) == 1 +*/ static SCM -scm_i_make_ratio (SCM numerator, SCM denominator) -#define FUNC_NAME "make-ratio" +scm_i_make_ratio_already_reduced (SCM numerator, SCM denominator) { - /* First make sure the arguments are proper. - */ - if (SCM_I_INUMP (denominator)) + /* Flip signs so that the denominator is positive. */ + if (scm_is_false (scm_positive_p (denominator))) { - if (scm_is_eq (denominator, SCM_INUM0)) + if (SCM_UNLIKELY (scm_is_eq (denominator, SCM_INUM0))) scm_num_overflow ("make-ratio"); - if (scm_is_eq (denominator, SCM_INUM1)) - return numerator; - } - else - { - if (!(SCM_BIGP(denominator))) - SCM_WRONG_TYPE_ARG (2, denominator); - } - if (!SCM_I_INUMP (numerator) && !SCM_BIGP (numerator)) - SCM_WRONG_TYPE_ARG (1, numerator); - - /* Then flip signs so that the denominator is positive. - */ - if (scm_is_true (scm_negative_p (denominator))) - { - numerator = scm_difference (numerator, SCM_UNDEFINED); - denominator = scm_difference (denominator, SCM_UNDEFINED); - } - - /* Now consider for each of the four fixnum/bignum combinations - whether the rational number is really an integer. - */ - if (SCM_I_INUMP (numerator)) - { - scm_t_inum x = SCM_I_INUM (numerator); - if (scm_is_eq (numerator, SCM_INUM0)) - return SCM_INUM0; - if (SCM_I_INUMP (denominator)) + else { - scm_t_inum y; - y = SCM_I_INUM (denominator); - if (x == y) - return SCM_INUM1; - if ((x % y) == 0) - return SCM_I_MAKINUM (x / y); + numerator = scm_difference (numerator, SCM_UNDEFINED); + denominator = scm_difference (denominator, SCM_UNDEFINED); } - else - { - /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative - of that value for the denominator, as a bignum. Apart from - that case, abs(bignum) > abs(inum) so inum/bignum is not an - integer. */ - if (x == SCM_MOST_NEGATIVE_FIXNUM - && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator), - - SCM_MOST_NEGATIVE_FIXNUM) == 0) - return SCM_I_MAKINUM(-1); - } } - else if (SCM_BIGP (numerator)) + + /* Check for the integer case */ + if (scm_is_eq (denominator, SCM_INUM1)) + return numerator; + + return scm_double_cell (scm_tc16_fraction, + SCM_UNPACK (numerator), + SCM_UNPACK (denominator), 0); +} + +static SCM scm_exact_integer_quotient (SCM x, SCM y); + +static SCM +scm_i_make_ratio (SCM numerator, SCM denominator) +#define FUNC_NAME "make-ratio" +{ + /* Make sure the arguments are proper */ + if (!SCM_LIKELY (SCM_I_INUMP (numerator) || SCM_BIGP (numerator))) + SCM_WRONG_TYPE_ARG (1, numerator); + else if (!SCM_LIKELY (SCM_I_INUMP (denominator) || SCM_BIGP (denominator))) + SCM_WRONG_TYPE_ARG (2, denominator); + else { - if (SCM_I_INUMP (denominator)) + SCM the_gcd = scm_gcd (numerator, denominator); + if (!(scm_is_eq (the_gcd, SCM_INUM1))) { - scm_t_inum yy = SCM_I_INUM (denominator); - if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy)) - return scm_divide (numerator, denominator); - } - else - { - if (scm_is_eq (numerator, denominator)) - return SCM_INUM1; - if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator), - SCM_I_BIG_MPZ (denominator))) - return scm_divide(numerator, denominator); + /* Reduce to lowest terms */ + numerator = scm_exact_integer_quotient (numerator, the_gcd); + denominator = scm_exact_integer_quotient (denominator, the_gcd); } + return scm_i_make_ratio_already_reduced (numerator, denominator); } - - /* No, it's a proper fraction. - */ - { - SCM divisor = scm_gcd (numerator, denominator); - if (!(scm_is_eq (divisor, SCM_INUM1))) - { - numerator = scm_divide (numerator, divisor); - denominator = scm_divide (denominator, divisor); - } - - return scm_double_cell (scm_tc16_fraction, - SCM_UNPACK (numerator), - SCM_UNPACK (denominator), 0); - } } #undef FUNC_NAME @@ -784,8 +746,9 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0, { if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x)))) return x; - return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED), - SCM_FRACTION_DENOMINATOR (x)); + return scm_i_make_ratio_already_reduced + (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED), + SCM_FRACTION_DENOMINATOR (x)); } else SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs); @@ -853,6 +816,83 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0, } #undef FUNC_NAME +/* scm_exact_integer_quotient returns 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). */ +static SCM +scm_exact_integer_quotient (SCM n, SCM d) +#define FUNC_NAME "exact-integer-quotient" +{ + if (SCM_LIKELY (SCM_I_INUMP (n))) + { + scm_t_inum nn = SCM_I_INUM (n); + if (SCM_LIKELY (SCM_I_INUMP (d))) + { + scm_t_inum dd = SCM_I_INUM (d); + if (SCM_UNLIKELY (dd == 0)) + scm_num_overflow ("exact-integer-quotient"); + else + { + scm_t_inum qq = nn / dd; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (qq); + } + } + else if (SCM_LIKELY (SCM_BIGP (d))) + { + /* n is an inum and d is a bignum. Given that d is known to + divide n evenly, there are only two possibilities: n is 0, + or else n is fixnum-min and d is abs(fixnum-min). */ + if (nn == 0) + return SCM_INUM0; + else + return SCM_I_MAKINUM (-1); + } + else + SCM_WRONG_TYPE_ARG (2, d); + } + else if (SCM_LIKELY (SCM_BIGP (n))) + { + if (SCM_LIKELY (SCM_I_INUMP (d))) + { + scm_t_inum dd = SCM_I_INUM (d); + if (SCM_UNLIKELY (dd == 0)) + scm_num_overflow ("exact-integer-quotient"); + else if (SCM_UNLIKELY (dd == 1)) + return n; + else + { + SCM q = scm_i_mkbig (); + if (dd > 0) + mpz_divexact_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), dd); + else + { + mpz_divexact_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), -dd); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + } + scm_remember_upto_here_1 (n); + return scm_i_normbig (q); + } + } + else if (SCM_LIKELY (SCM_BIGP (d))) + { + SCM q = scm_i_mkbig (); + mpz_divexact (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (n), + SCM_I_BIG_MPZ (d)); + scm_remember_upto_here_2 (n, d); + return scm_i_normbig (q); + } + else + SCM_WRONG_TYPE_ARG (2, d); + } + else + SCM_WRONG_TYPE_ARG (1, n); +} +#undef FUNC_NAME + /* two_valued_wta_dispatch_2 is a version of SCM_WTA_DISPATCH_2 for two-valued functions. It is called from primitive generics that take two arguments and return two values, when the core procedure is @@ -4636,6 +4676,23 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0, else /* return NaN for (0 ^ k) for negative k per R6RS */ return scm_nan (); } + /* Optimize the fraction case, to avoid needless reduction of + intermediate products to lowest terms. There will never be any + common factor to cancel out, so it would be a waste. */ + else if (SCM_FRACTIONP (n)) + { + if (scm_is_true (scm_positive_p (k))) + return scm_i_make_ratio_already_reduced + (scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k), + scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k)); + else + { + k = scm_difference (k, SCM_UNDEFINED); + return scm_i_make_ratio_already_reduced + (scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k), + scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k)); + } + } if (SCM_I_INUMP (k)) i2 = SCM_I_INUM (k); @@ -7282,8 +7339,9 @@ scm_difference (SCM x, SCM y) return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x), -SCM_COMPLEX_IMAG (x)); else if (SCM_FRACTIONP (x)) - return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED), - SCM_FRACTION_DENOMINATOR (x)); + return scm_i_make_ratio_already_reduced + (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED), + SCM_FRACTION_DENOMINATOR (x)); else SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference); } @@ -7825,14 +7883,14 @@ do_divide (SCM x, SCM y, int inexact) { if (inexact) return scm_from_double (1.0 / (double) xx); - else return scm_i_make_ratio (SCM_INUM1, x); + else return scm_i_make_ratio_already_reduced (SCM_INUM1, x); } } else if (SCM_BIGP (x)) { if (inexact) return scm_from_double (1.0 / scm_i_big2dbl (x)); - else return scm_i_make_ratio (SCM_INUM1, x); + else return scm_i_make_ratio_already_reduced (SCM_INUM1, x); } else if (SCM_REALP (x)) { @@ -7862,8 +7920,8 @@ do_divide (SCM x, SCM y, int inexact) } } else if (SCM_FRACTIONP (x)) - return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x), - SCM_FRACTION_NUMERATOR (x)); + return scm_i_make_ratio_already_reduced (SCM_FRACTION_DENOMINATOR (x), + SCM_FRACTION_NUMERATOR (x)); else SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide); } @@ -7960,32 +8018,9 @@ do_divide (SCM x, SCM y, int inexact) 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 - { - if (inexact) - return scm_from_double (scm_i_big2dbl (x) / (double) yy); - else return scm_i_make_ratio (x, y); - } + if (inexact) + return scm_from_double (scm_i_big2dbl (x) / (double) yy); + else return scm_i_make_ratio (x, y); } } else if (SCM_BIGP (y)) @@ -8003,21 +8038,7 @@ do_divide (SCM x, SCM y, int inexact) return scm_from_double (mpq_get_d (q)); } else - { - 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_i_make_ratio (x, y); } else if (SCM_REALP (y)) { @@ -8826,8 +8847,9 @@ SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 0, { if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z)))) return z; - return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z), SCM_UNDEFINED), - SCM_FRACTION_DENOMINATOR (z)); + return scm_i_make_ratio_already_reduced + (scm_difference (SCM_FRACTION_NUMERATOR (z), SCM_UNDEFINED), + SCM_FRACTION_DENOMINATOR (z)); } else SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude); @@ -8927,8 +8949,9 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0, mpq_init (frac); mpq_set_d (frac, val); - q = scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac)), - scm_i_mpz2num (mpq_denref (frac))); + q = scm_i_make_ratio_already_reduced + (scm_i_mpz2num (mpq_numref (frac)), + scm_i_mpz2num (mpq_denref (frac))); /* When scm_i_make_ratio throws, we leak the memory allocated for frac... diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index b1f3d8b..7cb4694 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4034,6 +4034,9 @@ (pass-if (eqv? -0.125 (expt -2 -3.0))) (pass-if (eqv? -0.125 (expt -2.0 -3.0))) (pass-if (eqv? 0.25 (expt 2.0 -2.0))) + (pass-if (eqv? 32/243 (expt 2/3 5))) + (pass-if (eqv? 243/32 (expt 2/3 -5))) + (pass-if (eqv? 32 (expt 1/2 -5))) (pass-if (test-eqv? (* -1.0+0.0i 12398 12398) (expt +12398i 2.0))) (pass-if (eqv-loosely? +i (expt -1 0.5))) (pass-if (eqv-loosely? +i (expt -1 1/2))) @@ -4304,6 +4307,9 @@ (pass-if (eqv? -1/8 (integer-expt -2 -3))) (pass-if (eqv? -0.125 (integer-expt -2.0 -3))) (pass-if (eqv? 0.25 (integer-expt 2.0 -2))) + (pass-if (eqv? 32/243 (integer-expt 2/3 5))) + (pass-if (eqv? 243/32 (integer-expt 2/3 -5))) + (pass-if (eqv? 32 (integer-expt 1/2 -5))) (pass-if (test-eqv? (* -1.0+0.0i 12398 12398) (integer-expt +12398.0i 2)))) -- 1.7.2.5