From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.org!not-for-mail From: Mark H Weaver Newsgroups: gmane.lisp.guile.devel Subject: Re: [PATCHES] Numeric improvements Date: Wed, 06 Mar 2013 19:16:05 -0500 Message-ID: <87fw089i9m.fsf@tines.lan> References: <87vc96ds5w.fsf@tines.lan> <87mwugd3o8.fsf@gnu.org> NNTP-Posting-Host: plane.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: ger.gmane.org 1362615421 9997 80.91.229.3 (7 Mar 2013 00:17:01 GMT) X-Complaints-To: usenet@ger.gmane.org NNTP-Posting-Date: Thu, 7 Mar 2013 00:17:01 +0000 (UTC) Cc: guile-devel@gnu.org To: ludo@gnu.org (Ludovic =?utf-8?Q?Court=C3=A8s?=) Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Thu Mar 07 01:17:21 2013 Return-path: Envelope-to: guile-devel@m.gmane.org Original-Received: from lists.gnu.org ([208.118.235.17]) by plane.gmane.org with esmtp (Exim 4.69) (envelope-from ) id 1UDOWK-0001AE-A2 for guile-devel@m.gmane.org; Thu, 07 Mar 2013 01:17:16 +0100 Original-Received: from localhost ([::1]:50303 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1UDOVy-00086W-Je for guile-devel@m.gmane.org; Wed, 06 Mar 2013 19:16:54 -0500 Original-Received: from eggs.gnu.org ([208.118.235.92]:55096) by lists.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1UDOVp-00085R-BW for guile-devel@gnu.org; Wed, 06 Mar 2013 19:16:49 -0500 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1UDOVj-0004O8-Ap for guile-devel@gnu.org; Wed, 06 Mar 2013 19:16:45 -0500 Original-Received: from world.peace.net ([96.39.62.75]:53060) by eggs.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1UDOVi-0004Ny-W0; Wed, 06 Mar 2013 19:16:39 -0500 Original-Received: from 209-6-91-212.c3-0.smr-ubr1.sbo-smr.ma.cable.rcn.com ([209.6.91.212] helo=tines.lan) by world.peace.net with esmtpsa (TLS1.0:DHE_RSA_AES_128_CBC_SHA1:16) (Exim 4.72) (envelope-from ) id 1UDOVJ-00085I-Aw; Wed, 06 Mar 2013 19:16:14 -0500 In-Reply-To: <87mwugd3o8.fsf@gnu.org> ("Ludovic \=\?utf-8\?Q\?Court\=C3\=A8s\=22'\?\= \=\?utf-8\?Q\?s\?\= message of "Wed, 06 Mar 2013 15:05:27 +0100") User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/24.3 (gnu/linux) X-detected-operating-system: by eggs.gnu.org: GNU/Linux 2.6.x X-Received-From: 96.39.62.75 X-BeenThere: guile-devel@gnu.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Developers list for Guile, the GNU extensibility library" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Original-Sender: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Xref: news.gmane.org gmane.lisp.guile.devel:15879 Archived-At: --=-=-= Content-Type: text/plain I pushed two of the patches: the 'gcd' patch and the one that verifies that FLT_RADIX is 2. Of the remaining 8 patches, I've attached improved versions of the first 4. Reviews welcome. Regards, Mark --=-=-= Content-Type: text/x-diff Content-Disposition: inline; filename=0001-Optimize-and-simplify-fractions-code.patch Content-Description: [PATCH 1/4] Optimize and simplify fractions code >From f32e8c5ffd789a6dbee48be74f5bbf32978382c3 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 3 Mar 2013 04:34:50 -0500 Subject: [PATCH 1/4] Optimize and simplify fractions code. * libguile/numbers.c (scm_exact_integer_quotient, scm_i_make_ratio_already_reduced): New static functions. (scm_i_make_ratio): Rewrite in terms of 'scm_i_make_ratio_already_reduced'. (scm_integer_expt): Optimize fraction case. (scm_abs, scm_magnitude, scm_difference, do_divide): Use 'scm_i_make_ratio_already_reduced'. * test-suite/tests/numbers.test (expt, integer-expt): Add tests. --- libguile/numbers.c | 244 ++++++++++++++++++++++++++--------------- test-suite/tests/numbers.test | 6 + 2 files changed, 159 insertions(+), 91 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index 393cf64..e63c17a 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -442,96 +442,56 @@ 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); +/* Make the ratio NUMERATOR/DENOMINATOR, where: + 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); + +/* Make the ratio NUMERATOR/DENOMINATOR */ +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 @@ -823,8 +783,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); @@ -892,6 +853,84 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0, } #undef FUNC_NAME +/* 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. */ +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 @@ -4675,6 +4714,26 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0, else /* return NaN for (0 ^ k) for negative k per R6RS */ return scm_nan (); } + else if (SCM_FRACTIONP (n)) + { + /* Optimize the fraction case by (a/b)^k ==> (a^k)/(b^k), to avoid + needless reduction of intermediate products to lowest terms. + If a and b have no common factors, then a^k and b^k have no + 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. */ + 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); @@ -7330,8 +7389,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); } @@ -7879,14 +7939,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)) { @@ -7916,8 +7976,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); } @@ -8880,8 +8940,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); @@ -8982,8 +9043,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 66aa01a..20b7eb1 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4044,6 +4044,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))) @@ -4317,6 +4320,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.10.4 --=-=-= Content-Type: text/x-diff Content-Disposition: inline; filename=0002-Add-round-ash-a-rounding-arithmetic-shift-operator.patch Content-Description: [PATCH 2/4] Add 'round-ash', a rounding arithmetic shift operator >From da13ef66508cd76b1c173ca355856b1e1b7667f0 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 3 Mar 2013 04:35:09 -0500 Subject: [PATCH 2/4] Add 'round-ash', a rounding arithmetic shift operator * libguile/numbers.c (left_shift_exact_integer, floor_right_shift_exact_integer, round_right_shift_exact_integer): New static functions. (scm_round_ash): New procedure. (scm_ash): Reimplement in terms of 'left_shift_exact_integer' and 'floor_right_shift_exact_integer'. * libguile/numbers.h: Add prototype for scm_round_ash. Rename the second argument of 'scm_ash' from 'cnt' to 'count'. * test-suite/tests/numbers.test (round-ash, ash): Add new unified testing framework for 'ash' and 'round-ash'. Previously, the tests for 'ash' were not very comprehensive; for example, they did not include a single test where the number to be shifted was a bignum. * doc/ref/api-data.texi (Bitwise Operations): Add documentation for 'round-ash'. Improve documentation for `ash'. --- doc/ref/api-data.texi | 44 +++++--- libguile/numbers.c | 226 +++++++++++++++++++++++++++-------------- libguile/numbers.h | 3 +- test-suite/tests/numbers.test | 114 +++++++++------------ 4 files changed, 235 insertions(+), 152 deletions(-) diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi index 9da17d8..b33270c 100644 --- a/doc/ref/api-data.texi +++ b/doc/ref/api-data.texi @@ -1686,19 +1686,16 @@ starts from 0 for the least significant bit. @end lisp @end deffn -@deffn {Scheme Procedure} ash n cnt -@deffnx {C Function} scm_ash (n, cnt) -Return @var{n} shifted left by @var{cnt} bits, or shifted right if -@var{cnt} is negative. This is an ``arithmetic'' shift. - -This is effectively a multiplication by @m{2^{cnt}, 2^@var{cnt}}, and -when @var{cnt} is negative it's a division, rounded towards negative -infinity. (Note that this is not the same rounding as @code{quotient} -does.) +@deffn {Scheme Procedure} ash n count +@deffnx {C Function} scm_ash (n, count) +Return @math{floor(@var{n} * 2^@var{count})}, but computed +more efficiently. @var{n} and @var{count} must be exact +integers. -With @var{n} viewed as an infinite precision twos complement, -@code{ash} means a left shift introducing zero bits, or a right shift -dropping bits. +With @var{n} viewed as an infinite-precision twos-complement +integer, @code{ash} means a left shift introducing zero bits +when @var{count} is positive, or a right shift dropping bits +when @var{count} is negative. This is an ``arithmetic'' shift. @lisp (number->string (ash #b1 3) 2) @result{} "1000" @@ -1709,6 +1706,29 @@ dropping bits. @end lisp @end deffn +@deffn {Scheme Procedure} round-ash n count +@deffnx {C Function} scm_round_ash (n, count) +Return @math{round(@var{n} * 2^@var{count})}, but computed +more efficiently. @var{n} and @var{count} must be exact +integers. + +With @var{n} viewed as an infinite-precision twos-complement +integer, @code{round-ash} means a left shift introducing zero +bits when @var{count} is positive, or a right shift rounding +to the nearest integer (with ties going to the nearest even +integer) when @var{count} is negative. This is a rounded +``arithmetic'' shift. + +@lisp +(number->string (round-ash #b1 3) 2) @result{} \"1000\" +(number->string (round-ash #b1010 -1) 2) @result{} \"101\" +(number->string (round-ash #b1010 -2) 2) @result{} \"10\" +(number->string (round-ash #b1011 -2) 2) @result{} \"11\" +(number->string (round-ash #b1101 -2) 2) @result{} \"11\" +(number->string (round-ash #b1110 -2) 2) @result{} \"100\" +@end lisp +@end deffn + @deffn {Scheme Procedure} logcount n @deffnx {C Function} scm_logcount (n) Return the number of bits in integer @var{n}. If @var{n} is diff --git a/libguile/numbers.c b/libguile/numbers.c index e63c17a..b99a04b 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -4791,19 +4791,119 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0, } #undef FUNC_NAME +/* Efficiently compute (N * 2^COUNT), + where N is an exact integer, and COUNT > 0. */ +static SCM +left_shift_exact_integer (SCM n, long count) +{ + if (SCM_I_INUMP (n)) + { + scm_t_inum nn = SCM_I_INUM (n); + + /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will always + overflow a non-zero fixnum. For smaller shifts we check the + bits going into positions above SCM_I_FIXNUM_BIT-1. If they're + all 0s for nn>=0, or all 1s for nn<0 then there's no overflow. + Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - count)". */ + + if (nn == 0) + return n; + else if (count < SCM_I_FIXNUM_BIT-1 && + ((scm_t_bits) (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - count)) + 1) + <= 1)) + return SCM_I_MAKINUM (nn << count); + else + { + SCM result = scm_i_inum2big (nn); + mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), + count); + return result; + } + } + else if (SCM_BIGP (n)) + { + SCM result = scm_i_mkbig (); + mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n), count); + scm_remember_upto_here_1 (n); + return result; + } + else + scm_syserror ("left_shift_exact_integer"); +} + +/* Efficiently compute floor (N / 2^COUNT), + where N is an exact integer and COUNT > 0. */ +static SCM +floor_right_shift_exact_integer (SCM n, long count) +{ + if (SCM_I_INUMP (n)) + { + scm_t_inum nn = SCM_I_INUM (n); + + if (count >= SCM_I_FIXNUM_BIT) + return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM (-1)); + else + return SCM_I_MAKINUM (SCM_SRS (nn, count)); + } + else if (SCM_BIGP (n)) + { + SCM result = scm_i_mkbig (); + mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n), + count); + scm_remember_upto_here_1 (n); + return scm_i_normbig (result); + } + else + scm_syserror ("floor_right_shift_exact_integer"); +} + +/* Efficiently compute round (N / 2^COUNT), + where N is an exact integer and COUNT > 0. */ +static SCM +round_right_shift_exact_integer (SCM n, long count) +{ + if (SCM_I_INUMP (n)) + { + if (count >= SCM_I_FIXNUM_BIT) + return SCM_INUM0; + else + { + scm_t_inum nn = SCM_I_INUM (n); + scm_t_inum qq = SCM_SRS (nn, count); + + if (0 == (nn & (1L << (count-1)))) + return SCM_I_MAKINUM (qq); /* round down */ + else if (nn & ((1L << (count-1)) - 1)) + return SCM_I_MAKINUM (qq + 1); /* round up */ + else + return SCM_I_MAKINUM ((~1L) & (qq + 1)); /* round to even */ + } + } + else if (SCM_BIGP (n)) + { + SCM q = scm_i_mkbig (); + + mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), count); + if (mpz_tstbit (SCM_I_BIG_MPZ (n), count-1) + && (mpz_odd_p (SCM_I_BIG_MPZ (q)) + || (mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1))) + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); + scm_remember_upto_here_1 (n); + return scm_i_normbig (q); + } + else + scm_syserror ("round_right_shift_exact_integer"); +} + SCM_DEFINE (scm_ash, "ash", 2, 0, 0, - (SCM n, SCM cnt), - "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n" - "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n" + (SCM n, SCM count), + "Return @math{floor(@var{n} * 2^@var{count})}.\n" + "@var{n} and @var{count} must be exact integers.\n" "\n" - "This is effectively a multiplication by 2^@var{cnt}, and when\n" - "@var{cnt} is negative it's a division, rounded towards negative\n" - "infinity. (Note that this is not the same rounding as\n" - "@code{quotient} does.)\n" - "\n" - "With @var{n} viewed as an infinite precision twos complement,\n" - "@code{ash} means a left shift introducing zero bits, or a right\n" - "shift dropping bits.\n" + "With @var{n} viewed as an infinite-precision twos-complement\n" + "integer, @code{ash} means a left shift introducing zero bits\n" + "when @var{count} is positive, or a right shift dropping bits\n" + "when @var{count} is negative. This is an ``arithmetic'' shift.\n" "\n" "@lisp\n" "(number->string (ash #b1 3) 2) @result{} \"1000\"\n" @@ -4814,79 +4914,57 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_ash { - long bits_to_shift; - bits_to_shift = scm_to_long (cnt); - - if (SCM_I_INUMP (n)) + if (SCM_I_INUMP (n) || SCM_BIGP (n)) { - scm_t_inum nn = SCM_I_INUM (n); + long bits_to_shift = scm_to_long (count); if (bits_to_shift > 0) - { - /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always - overflow a non-zero fixnum. For smaller shifts we check the - bits going into positions above SCM_I_FIXNUM_BIT-1. If they're - all 0s for nn>=0, or all 1s for nn<0 then there's no overflow. - Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - - bits_to_shift)". */ - - if (nn == 0) - return n; - - if (bits_to_shift < SCM_I_FIXNUM_BIT-1 - && ((scm_t_bits) - (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1) - <= 1)) - { - return SCM_I_MAKINUM (nn << bits_to_shift); - } - else - { - SCM result = scm_i_inum2big (nn); - mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), - bits_to_shift); - return result; - } - } + return left_shift_exact_integer (n, bits_to_shift); + else if (SCM_LIKELY (bits_to_shift < 0)) + return floor_right_shift_exact_integer (n, -bits_to_shift); else - { - bits_to_shift = -bits_to_shift; - if (bits_to_shift >= SCM_LONG_BIT) - return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM(-1)); - else - return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift)); - } - + return n; } - else if (SCM_BIGP (n)) - { - SCM result; + else + SCM_WRONG_TYPE_ARG (SCM_ARG1, n); +} +#undef FUNC_NAME - if (bits_to_shift == 0) - return n; +SCM_DEFINE (scm_round_ash, "round-ash", 2, 0, 0, + (SCM n, SCM count), + "Return @math{round(@var{n} * 2^@var{count})}.\n" + "@var{n} and @var{count} must be exact integers.\n" + "\n" + "With @var{n} viewed as an infinite-precision twos-complement\n" + "integer, @code{round-ash} means a left shift introducing zero\n" + "bits when @var{count} is positive, or a right shift rounding\n" + "to the nearest integer (with ties going to the nearest even\n" + "integer) when @var{count} is negative. This is a rounded\n" + "``arithmetic'' shift.\n" + "\n" + "@lisp\n" + "(number->string (round-ash #b1 3) 2) @result{} \"1000\"\n" + "(number->string (round-ash #b1010 -1) 2) @result{} \"101\"\n" + "(number->string (round-ash #b1010 -2) 2) @result{} \"10\"\n" + "(number->string (round-ash #b1011 -2) 2) @result{} \"11\"\n" + "(number->string (round-ash #b1101 -2) 2) @result{} \"11\"\n" + "(number->string (round-ash #b1110 -2) 2) @result{} \"100\"\n" + "@end lisp") +#define FUNC_NAME s_scm_round_ash +{ + if (SCM_I_INUMP (n) || SCM_BIGP (n)) + { + long bits_to_shift = scm_to_long (count); - result = scm_i_mkbig (); - if (bits_to_shift >= 0) - { - mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n), - bits_to_shift); - return result; - } + if (bits_to_shift > 0) + return left_shift_exact_integer (n, bits_to_shift); + else if (SCM_LIKELY (bits_to_shift < 0)) + return round_right_shift_exact_integer (n, -bits_to_shift); else - { - /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so - we have to allocate a bignum even if the result is going to be a - fixnum. */ - mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n), - -bits_to_shift); - return scm_i_normbig (result); - } - + return n; } else - { - SCM_WRONG_TYPE_ARG (SCM_ARG1, n); - } + SCM_WRONG_TYPE_ARG (SCM_ARG1, n); } #undef FUNC_NAME diff --git a/libguile/numbers.h b/libguile/numbers.h index 2c8b260..912f287 100644 --- a/libguile/numbers.h +++ b/libguile/numbers.h @@ -206,7 +206,8 @@ SCM_API SCM scm_logbit_p (SCM n1, SCM n2); SCM_API SCM scm_lognot (SCM n); SCM_API SCM scm_modulo_expt (SCM n, SCM k, SCM m); SCM_API SCM scm_integer_expt (SCM z1, SCM z2); -SCM_API SCM scm_ash (SCM n, SCM cnt); +SCM_API SCM scm_ash (SCM n, SCM count); +SCM_API SCM scm_round_ash (SCM n, SCM count); SCM_API SCM scm_bit_extract (SCM n, SCM start, SCM end); SCM_API SCM scm_logcount (SCM n); SCM_API SCM scm_integer_length (SCM n); diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index 20b7eb1..8dab29c 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -201,71 +201,6 @@ (eqv? -2305843009213693953 (1- -2305843009213693952)))) ;;; -;;; ash -;;; - -(with-test-prefix "ash" - - (pass-if "documented?" - (documented? ash)) - - (pass-if (eqv? 0 (ash 0 0))) - (pass-if (eqv? 0 (ash 0 1))) - (pass-if (eqv? 0 (ash 0 1000))) - (pass-if (eqv? 0 (ash 0 -1))) - (pass-if (eqv? 0 (ash 0 -1000))) - - (pass-if (eqv? 1 (ash 1 0))) - (pass-if (eqv? 2 (ash 1 1))) - (pass-if (eqv? 340282366920938463463374607431768211456 (ash 1 128))) - (pass-if (eqv? 0 (ash 1 -1))) - (pass-if (eqv? 0 (ash 1 -1000))) - - (pass-if (eqv? -1 (ash -1 0))) - (pass-if (eqv? -2 (ash -1 1))) - (pass-if (eqv? -340282366920938463463374607431768211456 (ash -1 128))) - (pass-if (eqv? -1 (ash -1 -1))) - (pass-if (eqv? -1 (ash -1 -1000))) - - (pass-if (eqv? -3 (ash -3 0))) - (pass-if (eqv? -6 (ash -3 1))) - (pass-if (eqv? -1020847100762815390390123822295304634368 (ash -3 128))) - (pass-if (eqv? -2 (ash -3 -1))) - (pass-if (eqv? -1 (ash -3 -1000))) - - (pass-if (eqv? -6 (ash -23 -2))) - - (pass-if (eqv? most-positive-fixnum (ash most-positive-fixnum 0))) - (pass-if (eqv? (* 2 most-positive-fixnum) (ash most-positive-fixnum 1))) - (pass-if (eqv? (* 4 most-positive-fixnum) (ash most-positive-fixnum 2))) - (pass-if - (eqv? (* most-positive-fixnum 340282366920938463463374607431768211456) - (ash most-positive-fixnum 128))) - (pass-if (eqv? (quotient most-positive-fixnum 2) - (ash most-positive-fixnum -1))) - (pass-if (eqv? 0 (ash most-positive-fixnum -1000))) - - (let ((mpf4 (quotient most-positive-fixnum 4))) - (pass-if (eqv? (* 2 mpf4) (ash mpf4 1))) - (pass-if (eqv? (* 4 mpf4) (ash mpf4 2))) - (pass-if (eqv? (* 8 mpf4) (ash mpf4 3)))) - - (pass-if (eqv? most-negative-fixnum (ash most-negative-fixnum 0))) - (pass-if (eqv? (* 2 most-negative-fixnum) (ash most-negative-fixnum 1))) - (pass-if (eqv? (* 4 most-negative-fixnum) (ash most-negative-fixnum 2))) - (pass-if - (eqv? (* most-negative-fixnum 340282366920938463463374607431768211456) - (ash most-negative-fixnum 128))) - (pass-if (eqv? (quotient-floor most-negative-fixnum 2) - (ash most-negative-fixnum -1))) - (pass-if (eqv? -1 (ash most-negative-fixnum -1000))) - - (let ((mnf4 (quotient-floor most-negative-fixnum 4))) - (pass-if (eqv? (* 2 mnf4) (ash mnf4 1))) - (pass-if (eqv? (* 4 mnf4) (ash mnf4 2))) - (pass-if (eqv? (* 8 mnf4) (ash mnf4 3))))) - -;;; ;;; exact? ;;; @@ -4904,3 +4839,52 @@ round-quotient round-remainder valid-round-answer?))) + +;;; +;;; ash +;;; round-ash +;;; + +(let () + (define (test-ash-variant name ash-variant round-variant) + (with-test-prefix name + (define (test n count) + (pass-if (list n count) + (eqv? (ash-variant n count) + (round-variant (* n (expt 2 count)))))) + + (pass-if "documented?" + (documented? ash-variant)) + + (for-each (lambda (n) + (for-each (lambda (count) (test n count)) + '(-1000 -3 -2 -1 0 1 2 3 1000))) + (list 0 1 3 23 -1 -3 -23 + fixnum-max + (1+ fixnum-max) + (1- fixnum-max) + (* fixnum-max 4) + (quotient fixnum-max 4) + fixnum-min + (1+ fixnum-min) + (1- fixnum-min) + (* fixnum-min 4) + (quotient fixnum-min 4))) + + (do ((count -2 (1- count)) + (vals '(1 3 5 7 9 11) + (map (lambda (n) (* 2 n)) vals))) + ((> (car vals) (* 2 fixnum-max)) 'done) + (for-each (lambda (n) + (test n count) + (test (- n) count)) + vals)) + + ;; Test rounding + (for-each (lambda (base) + (for-each (lambda (offset) (test (+ base offset) -3)) + '(#b11001 #b11100 #b11101 #b10001 #b10100 #b10101))) + (list 0 64 -64 (* 64 fixnum-max) (* 64 fixnum-min))))) + + (test-ash-variant 'ash ash floor) + (test-ash-variant 'round-ash round-ash round)) -- 1.7.10.4 --=-=-= Content-Type: text/x-diff Content-Disposition: inline; filename=0003-Simplify-and-improve-scm_i_big2dbl-and-add-scm_i_big.patch Content-Description: [PATCH 3/4] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp >From f48e23eba996566adbe112731ceda4761070fbce Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 3 Mar 2013 04:58:55 -0500 Subject: [PATCH 3/4] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp * libguile/numbers.c (scm_i_big2dbl_2exp): New static function. (scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp', with proper rounding. * test-suite/tests/numbers.test ("inexact->exact"): Add tests. --- libguile/numbers.c | 101 +++++++++++++++-------------------------- test-suite/tests/numbers.test | 45 ++++++++++++++++++ 2 files changed, 81 insertions(+), 65 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index b99a04b..49b4a50 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -330,81 +330,52 @@ scm_i_dbl2num (double u) return scm_i_dbl2big (u); } -/* scm_i_big2dbl() rounds to the closest representable double, in accordance - with R5RS exact->inexact. +static SCM round_right_shift_exact_integer (SCM n, long count); - The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits - (ie. truncate towards zero), then adjust to get the closest double by - examining the next lower bit and adding 1 (to the absolute value) if - necessary. +/* scm_i_big2dbl_2exp() is like frexp for bignums: it converts the + bignum b into a normalized significand and exponent such that + b = significand * 2^exponent and 1/2 <= abs(significand) < 1. + The return value is the significand rounded to the closest + representable double, and the exponent is placed into *expon_p. + If b is zero, then the returned exponent and significand are both + zero. */ - Bignums exactly half way between representable doubles are rounded to the - next higher absolute value (ie. away from zero). This seems like an - adequate interpretation of R5RS "numerically closest", and it's easier - and faster than a full "nearest-even" style. - - The bit test must be done on the absolute value of the mpz_t, which means - we need to use mpz_getlimbn. mpz_tstbit is not right, it treats - negatives as twos complement. - - In GMP before 4.2, mpz_get_d rounding was unspecified. It ended up - following the hardware rounding mode, but applied to the absolute - value of the mpz_t operand. This is not what we want so we put the - high DBL_MANT_DIG bits into a temporary. Starting with GMP 4.2 - (released in March 2006) mpz_get_d now always truncates towards zero. - - ENHANCE-ME: The temporary init+clear to force the rounding in GMP - before 4.2 is a slowdown. It'd be faster to pick out the relevant - high bits with mpz_getlimbn. */ - -double -scm_i_big2dbl (SCM b) +static double +scm_i_big2dbl_2exp (SCM b, long *expon_p) { - double result; - size_t bits; - - bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2); - -#if 1 - { - /* For GMP earlier than 4.2, force truncation towards zero */ - - /* FIXME: DBL_MANT_DIG is the number of base-`FLT_RADIX' digits, - _not_ the number of bits, so this code will break badly on a - system with non-binary doubles. */ - - mpz_t tmp; - if (bits > DBL_MANT_DIG) - { - size_t shift = bits - DBL_MANT_DIG; - mpz_init2 (tmp, DBL_MANT_DIG); - mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift); - result = ldexp (mpz_get_d (tmp), shift); - mpz_clear (tmp); - } - else - { - result = mpz_get_d (SCM_I_BIG_MPZ (b)); - } - } -#else - /* GMP 4.2 or later */ - result = mpz_get_d (SCM_I_BIG_MPZ (b)); -#endif + size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2); + size_t shift = 0; if (bits > DBL_MANT_DIG) { - unsigned long pos = bits - DBL_MANT_DIG - 1; - /* test bit number "pos" in absolute value */ - if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS) - & ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS))) + shift = bits - DBL_MANT_DIG; + b = round_right_shift_exact_integer (b, shift); + if (SCM_I_INUMP (b)) { - result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1); + int expon; + double signif = frexp (SCM_I_INUM (b), &expon); + *expon_p = expon + shift; + return signif; } } - scm_remember_upto_here_1 (b); - return result; + { + long expon; + double signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b)); + scm_remember_upto_here_1 (b); + *expon_p = expon + shift; + return signif; + } +} + +/* scm_i_big2dbl() rounds to the closest representable double, + in accordance with R5RS exact->inexact. */ +double +scm_i_big2dbl (SCM b) +{ + long expon; + double signif = scm_i_big2dbl_2exp (b, &expon); + return ldexp (signif, expon); } SCM diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index 8dab29c..e3a1099 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4220,6 +4220,51 @@ (pass-if-exception "nan" exception:out-of-range (inexact->exact +nan.0)) + + (pass-if-equal "round up to odd" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111000101 -> + ;; 11111111111111111111111111111111111111111111111111001000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000101)))) + + (pass-if-equal "round down to odd" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111001011 -> + ;; 11111111111111111111111111111111111111111111111111001000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001011)))) + + (pass-if-equal "round tie up to even" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111011100 -> + ;; 11111111111111111111111111111111111111111111111111100000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b100000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b011100)))) + + (pass-if-equal "round tie down to even" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111000100 -> + ;; 11111111111111111111111111111111111111111111111111000000 + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000000) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000100)))) + + (pass-if-equal "round tie up to next power of two" + ;; ===================================================== + ;; 11111111111111111111111111111111111111111111111111111100 -> + ;; 100000000000000000000000000000000000000000000000000000000 + (expt 2 (+ dbl-mant-dig 3)) + (inexact->exact + (exact->inexact + (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b111100)))) (with-test-prefix "2.0**i to exact and back" (do ((i 0 (1+ i)) -- 1.7.10.4 --=-=-= Content-Type: text/x-diff Content-Disposition: inline; filename=0004-Optimize-logarithms-using-scm_i_big2dbl_2exp.patch Content-Description: [PATCH 4/4] Optimize logarithms using scm_i_big2dbl_2exp >From b21e3d03049e986b0ae6e0ad59aa79c6c38de842 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Sun, 3 Mar 2013 05:02:53 -0500 Subject: [PATCH 4/4] Optimize logarithms using scm_i_big2dbl_2exp * libguile/numbers.c (log_of_exact_integer_with_size): Removed. (log_of_exact_integer): Handle bignums too large to fit in a double using 'scm_i_big2dbl_2exp' instead of 'scm_integer_length' and 'scm_ash'. (log_of_fraction): Use 'log_of_exact_integer' instead of 'log_of_exact_integer_with_size'. --- libguile/numbers.c | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index 49b4a50..0b13d69 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -9576,26 +9576,20 @@ log_of_shifted_double (double x, long shift) return scm_c_make_rectangular (ans, M_PI); } -/* Returns log(n), for exact integer n of integer-length size */ -static SCM -log_of_exact_integer_with_size (SCM n, long size) -{ - long shift = size - 2 * scm_dblprec[0]; - - if (shift > 0) - return log_of_shifted_double - (scm_to_double (scm_ash (n, scm_from_long(-shift))), - shift); - else - return log_of_shifted_double (scm_to_double (n), 0); -} - /* Returns log(n), for exact integer n */ static SCM log_of_exact_integer (SCM n) { - return log_of_exact_integer_with_size - (n, scm_to_long (scm_integer_length (n))); + if (SCM_I_INUMP (n)) + return log_of_shifted_double (SCM_I_INUM (n), 0); + else if (SCM_BIGP (n)) + { + long expon; + double signif = scm_i_big2dbl_2exp (n, &expon); + return log_of_shifted_double (signif, expon); + } + else + scm_wrong_type_arg ("log_of_exact_integer", SCM_ARG1, n); } /* Returns log(n/d), for exact non-zero integers n and d */ @@ -9606,8 +9600,8 @@ log_of_fraction (SCM n, SCM d) long d_size = scm_to_long (scm_integer_length (d)); if (abs (n_size - d_size) > 1) - return (scm_difference (log_of_exact_integer_with_size (n, n_size), - log_of_exact_integer_with_size (d, d_size))); + return (scm_difference (log_of_exact_integer (n), + log_of_exact_integer (d))); else if (scm_is_false (scm_negative_p (n))) return scm_from_double (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d)))); -- 1.7.10.4 --=-=-=--