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: [PATCH] Improved log/log10, add round-ash Date: Tue, 15 Feb 2011 11:12:02 -0500 Message-ID: <871v395lbh.fsf@netris.org> NNTP-Posting-Host: lo.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: dough.gmane.org 1297786386 28631 80.91.229.12 (15 Feb 2011 16:13:06 GMT) X-Complaints-To: usenet@dough.gmane.org NNTP-Posting-Date: Tue, 15 Feb 2011 16:13:06 +0000 (UTC) To: guile-devel@gnu.org Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Tue Feb 15 17:13:00 2011 Return-path: Envelope-to: guile-devel@m.gmane.org Original-Received: from lists.gnu.org ([199.232.76.165]) by lo.gmane.org with esmtp (Exim 4.69) (envelope-from ) id 1PpNWB-00020b-GF for guile-devel@m.gmane.org; Tue, 15 Feb 2011 17:12:59 +0100 Original-Received: from localhost ([127.0.0.1]:36531 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PpNWA-0006ey-EN for guile-devel@m.gmane.org; Tue, 15 Feb 2011 11:12:46 -0500 Original-Received: from [140.186.70.92] (port=39821 helo=eggs.gnu.org) by lists.gnu.org with esmtp (Exim 4.43) id 1PpNVu-0006Yn-GG for guile-devel@gnu.org; Tue, 15 Feb 2011 11:12:40 -0500 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1PpNVe-0000vV-AG for guile-devel@gnu.org; Tue, 15 Feb 2011 11:12:30 -0500 Original-Received: from world.peace.net ([216.204.32.208]:51273) by eggs.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1PpNVd-0000tu-DO for guile-devel@gnu.org; Tue, 15 Feb 2011 11:12:13 -0500 Original-Received: from ip68-9-118-38.ri.ri.cox.net ([68.9.118.38] helo=freedomincluded) by world.peace.net with esmtpa (Exim 4.69) (envelope-from ) id 1PpNVT-0007Wd-8w; Tue, 15 Feb 2011 11:12:03 -0500 Original-Received: from mhw by freedomincluded with local (Exim 4.69) (envelope-from ) id 1PpNVS-0003dQ-Hp; Tue, 15 Feb 2011 11:12:02 -0500 User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/23.2 (gnu/linux) X-detected-operating-system: by eggs.gnu.org: GNU/Linux 2.6 (newer, 3) X-Received-From: 216.204.32.208 X-BeenThere: guile-devel@gnu.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Developers list for Guile, the GNU extensibility library" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Original-Sender: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Errors-To: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Xref: news.gmane.org gmane.lisp.guile.devel:11649 Archived-At: --=-=-= The first patch here is a trivial comment fix for numbers.test. The second patch fixes some serious problems in the logarithm functions when applied to large integers, large or small fractions, or fractions with large numerators and denominators which are close to 1. I strongly recommend including it in 2.0.0. The third patch contains round-ash, the new rounding arithmetic shift operator. I think it's useful and would like to see it in 2.0, but will understand if you wish otherwise. This version has been rebased to the current trunk. Best, Mark --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0001-Fix-comment-above-number-theoretic-division-tests.patch Content-Description: Fix comment above number-theoretic division tests >From 4716bd257fb2d709c44d717ebb9ad2d8f4b50b8e Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Mon, 14 Feb 2011 18:18:52 -0500 Subject: [PATCH 1/3] Fix comment above number-theoretic division tests * test-suite/tests/numbers.test: Fix comment. --- test-suite/tests/numbers.test | 14 ++++++++++++++ 1 files changed, 14 insertions(+), 0 deletions(-) diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index 1f2ee03..9e9728f 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4512,12 +4512,26 @@ ;;; +;;; Tests for number-theoretic division operators: +;;; ;;; euclidean/ ;;; euclidean-quotient ;;; euclidean-remainder +;;; floor/ +;;; floor-quotient +;;; floor-remainder +;;; ceiling/ +;;; ceiling-quotient +;;; ceiling-remainder +;;; truncate/ +;;; truncate-quotient +;;; truncate-remainder ;;; centered/ ;;; centered-quotient ;;; centered-remainder +;;; round/ +;;; round-quotient +;;; round-remainder ;;; (with-test-prefix "Number-theoretic division" -- 1.7.1 --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0002-Improvements-to-log-and-log10.patch Content-Description: Improvements to `log' and `log10' >From 057df863b37bbc71010e647608a9daa7ae876afc Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Tue, 15 Feb 2011 10:37:03 -0500 Subject: [PATCH 2/3] Improvements to `log' and `log10' * libguile/numbers.c (log_of_shifted_double, log_of_exact_integer, log_of_exact_integer_with_size, log_of_fraction): New internal static functions used by scm_log and scm_log10. (scm_log, scm_log10): Robustly handle large integers, large and small fractions, and fractions close to 1. Previously, computing logarithms of fractions close to 1 yielded grossly inaccurate results, and the other cases yielded infinities even though the answer could easily fit in a double. (log -0.0) now returns -inf.0+i, where previously it returned -inf.0. (log 0) now throws a numerical overflow exception, where previously it returned -inf.0. (log 0.0) still returns -inf.0. Analogous changes made to `log10'. * test-suite/tests/numbers.test (log, log10): Add tests. --- libguile/numbers.c | 108 ++++++++++++++++++++++++++++++++++------- test-suite/tests/numbers.test | 80 +++++++++++++++++++++++------- 2 files changed, 152 insertions(+), 36 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index 7c4ea1b..d0aacb7 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -111,6 +111,7 @@ typedef scm_t_signed_bits scm_t_inum; static SCM flo0; static SCM exactly_one_half; +static SCM flo_log10e; #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0) @@ -9372,6 +9373,62 @@ scm_is_number (SCM z) } +/* Returns log(x * 2^shift) */ +static SCM +log_of_shifted_double (double x, long shift) +{ + double ans = log (fabs (x)) + shift * M_LN2; + + if (x > 0.0 || double_is_non_negative_zero (x)) + return scm_from_double (ans); + else + 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 of integer-length size */ +static SCM +log_of_exact_integer (SCM n) +{ + return log_of_exact_integer_with_size + (n, scm_to_long (scm_integer_length (n))); +} + +/* Returns log(n/d), for exact non-zero integers n and d */ +static SCM +log_of_fraction (SCM n, SCM d) +{ + long n_size = scm_to_long (scm_integer_length (n)); + 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))); + else if (scm_is_false (scm_negative_p (n))) + return scm_from_double + (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d)))); + else + return scm_c_make_rectangular + (log1p (scm_to_double (scm_divide2real + (scm_difference (scm_abs (n), d), + d))), + M_PI); +} + + /* In the following functions we dispatch to the real-arg funcs like log() when we know the arg is real, instead of just handing everything to clog() for instance. This is in case clog() doesn't optimize for a @@ -9394,17 +9451,21 @@ SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0, atan2 (im, re)); #endif } - else if (SCM_NUMBERP (z)) + else if (SCM_REALP (z)) + return log_of_shifted_double (SCM_REAL_VALUE (z), 0); + else if (SCM_I_INUMP (z)) { - /* ENHANCE-ME: When z is a bignum the logarithm will fit a double - although the value itself overflows. */ - double re = scm_to_double (z); - double l = log (fabs (re)); - if (re >= 0.0) - return scm_from_double (l); - else - return scm_c_make_rectangular (l, M_PI); +#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO + if (scm_is_eq (z, SCM_INUM0)) + scm_num_overflow (s_scm_log); +#endif + return log_of_shifted_double (SCM_I_INUM (z), 0); } + else if (SCM_BIGP (z)) + return log_of_exact_integer (z); + else if (SCM_FRACTIONP (z)) + return log_of_fraction (SCM_FRACTION_NUMERATOR (z), + SCM_FRACTION_DENOMINATOR (z)); else SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log); } @@ -9431,17 +9492,27 @@ SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0, M_LOG10E * atan2 (im, re)); #endif } - else if (SCM_NUMBERP (z)) + else if (SCM_REALP (z) || SCM_I_INUMP (z)) { - /* ENHANCE-ME: When z is a bignum the logarithm will fit a double - although the value itself overflows. */ - double re = scm_to_double (z); - double l = log10 (fabs (re)); - if (re >= 0.0) - return scm_from_double (l); - else - return scm_c_make_rectangular (l, M_LOG10E * M_PI); +#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO + if (scm_is_eq (z, SCM_INUM0)) + scm_num_overflow (s_scm_log10); +#endif + { + double re = scm_to_double (z); + double l = log10 (fabs (re)); + if (re > 0.0 || double_is_non_negative_zero (re)) + return scm_from_double (l); + else + return scm_c_make_rectangular (l, M_LOG10E * M_PI); + } } + else if (SCM_BIGP (z)) + return scm_product (flo_log10e, log_of_exact_integer (z)); + else if (SCM_FRACTIONP (z)) + return scm_product (flo_log10e, + log_of_fraction (SCM_FRACTION_NUMERATOR (z), + SCM_FRACTION_DENOMINATOR (z))); else SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10); } @@ -9536,6 +9607,7 @@ scm_init_numbers () scm_add_feature ("complex"); scm_add_feature ("inexact"); flo0 = scm_from_double (0.0); + flo_log10e = scm_from_double (M_LOG10E); /* determine floating point precision */ for (i=2; i <= SCM_MAX_DBL_RADIX; ++i) diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test index 9e9728f..cb582ed 100644 --- a/test-suite/tests/numbers.test +++ b/test-suite/tests/numbers.test @@ -4323,14 +4323,36 @@ (log)) (pass-if-exception "two args" exception:wrong-num-args (log 123 456)) - - (pass-if (negative-infinity? (log 0))) - (pass-if (negative-infinity? (log 0.0))) - (pass-if (eqv? 0.0 (log 1))) - (pass-if (eqv? 0.0 (log 1.0))) - (pass-if (eqv-loosely? 1.0 (log const-e))) - (pass-if (eqv-loosely? 2.0 (log const-e^2))) - (pass-if (eqv-loosely? -1.0 (log const-1/e))) + (pass-if-exception "(log 0)" exception:numerical-overflow + (log 0)) + + (pass-if (test-eqv? -inf.0 (log 0.0))) + (pass-if (test-eqv? +inf.0 (log +inf.0))) + (pass-if (test-eqv? -inf.0+3.14159265358979i (log -0.0))) + (pass-if (test-eqv? +inf.0+3.14159265358979i (log -inf.0))) + (pass-if (test-eqv? 0.0 (log 1 ))) + (pass-if (test-eqv? 0.0 (log 1.0))) + (pass-if (test-eqv? 1.0 (log const-e))) + (pass-if (test-eqv? 2.0 (log const-e^2))) + (pass-if (test-eqv? -1.0 (log const-1/e))) + (pass-if (test-eqv? -1.0+3.14159265358979i (log (- const-1/e)))) + (pass-if (test-eqv? 2.30258509299405 (log 10))) + (pass-if (test-eqv? 2.30258509299405+3.14159265358979i (log -10))) + + (pass-if (test-eqv? 1.0+0.0i (log (+ const-e +0.0i)))) + (pass-if (test-eqv? 1.0-0.0i (log (+ const-e -0.0i)))) + + (pass-if (eqv-loosely? 230258.509299405 (log (expt 10 100000)))) + (pass-if (eqv-loosely? -230258.509299405 (log (expt 10 -100000)))) + (pass-if (eqv-loosely? 230257.410687116 (log (/ (expt 10 100000) 3)))) + (pass-if (eqv-loosely? 230258.509299405+3.14159265358979i + (log (- (expt 10 100000))))) + (pass-if (eqv-loosely? -230258.509299405+3.14159265358979i + (log (- (expt 10 -100000))))) + (pass-if (eqv-loosely? 230257.410687116+3.14159265358979i + (log (- (/ (expt 10 100000) 3))))) + (pass-if (test-eqv? 3.05493636349961e-151 + (log (/ (1+ (expt 2 500)) (expt 2 500))))) (pass-if (eqv-loosely? 1.0+1.57079i (log 0+2.71828i))) (pass-if (eqv-loosely? 1.0-1.57079i (log 0-2.71828i))) @@ -4350,20 +4372,42 @@ (log10)) (pass-if-exception "two args" exception:wrong-num-args (log10 123 456)) - - (pass-if (negative-infinity? (log10 0))) - (pass-if (negative-infinity? (log10 0.0))) - (pass-if (eqv? 0.0 (log10 1))) - (pass-if (eqv? 0.0 (log10 1.0))) - (pass-if (eqv-loosely? 1.0 (log10 10.0))) - (pass-if (eqv-loosely? 2.0 (log10 100.0))) - (pass-if (eqv-loosely? -1.0 (log10 0.1))) + (pass-if-exception "(log10 0)" exception:numerical-overflow + (log10 0)) + + (pass-if (test-eqv? -inf.0 (log10 0.0))) + (pass-if (test-eqv? +inf.0 (log10 +inf.0))) + (pass-if (test-eqv? -inf.0+1.36437635384184i (log10 -0.0))) + (pass-if (test-eqv? +inf.0+1.36437635384184i (log10 -inf.0))) + (pass-if (test-eqv? 0.0 (log10 1 ))) + (pass-if (test-eqv? 0.0 (log10 1.0))) + (pass-if (test-eqv? 1.0 (log10 10 ))) + (pass-if (test-eqv? 1.0 (log10 10.0))) + (pass-if (test-eqv? 2.0 (log10 100.0))) + (pass-if (test-eqv? -1.0 (log10 0.1))) + (pass-if (test-eqv? -1.0+1.36437635384184i (log10 -0.1))) + (pass-if (test-eqv? 1.0+1.36437635384184i (log10 -10 ))) + + (pass-if (test-eqv? 1.0+0.0i (log10 10.0+0.0i))) + (pass-if (test-eqv? 1.0-0.0i (log10 10.0-0.0i))) + + (pass-if (eqv-loosely? 100000.0 (log10 (expt 10 100000)))) + (pass-if (eqv-loosely? -100000.0 (log10 (expt 10 -100000)))) + (pass-if (eqv-loosely? 99999.5228787453 (log10 (/ (expt 10 100000) 3)))) + (pass-if (eqv-loosely? 100000.0+1.36437635384184i + (log10 (- (expt 10 100000))))) + (pass-if (eqv-loosely? -100000.0+1.36437635384184i + (log10 (- (expt 10 -100000))))) + (pass-if (eqv-loosely? 99999.5228787453+1.36437635384184i + (log10 (- (/ (expt 10 100000) 3))))) + (pass-if (test-eqv? 1.32674200523347e-151 + (log10 (/ (1+ (expt 2 500)) (expt 2 500))))) (pass-if (eqv-loosely? 1.0+0.68218i (log10 0+10.0i))) (pass-if (eqv-loosely? 1.0-0.68218i (log10 0-10.0i))) - (pass-if (eqv-loosely? 0.0+1.36437i (log10 -1))) - (pass-if (eqv-loosely? 1.0+1.36437i (log10 -10))) + (pass-if (eqv-loosely? 0.0+1.36437i (log10 -1))) + (pass-if (eqv-loosely? 1.0+1.36437i (log10 -10))) (pass-if (eqv-loosely? 2.0+1.36437i (log10 -100)))) ;;; -- 1.7.1 --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0003-Add-round-ash-a-rounding-arithmetic-shift-operator.patch Content-Description: Add `round-ash', a rounding arithmetic shift operator >From 16ca1ed41ecb23a69fbd421f2f7d380bcdc82fd0 Mon Sep 17 00:00:00 2001 From: Mark H Weaver Date: Tue, 15 Feb 2011 07:34:54 -0500 Subject: [PATCH 3/3] 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 internal static functions to efficiently compute n * 2^count, floor(n / 2^count), and round(n / 2^count), respectively, for exact integer N and positive exact integer COUNT. Used to combine the implementations of `ash' and `round-ash' with minimal code duplication. (scm_ash): Reimplement in terms of left_shift_exact_integer and floor_right_shift_exact_integer. Note that this function efficiently computes floor(n * 2^count). (scm_round_ash): New procedure to efficiently compute round(n * 2^count), for exact integers N and COUNT. Implemented in terms of left_shift_exact_integer and round_right_shift_exact_integer. * libguile/numbers.h: Add prototype for scm_round_ash. Change the name of scm_ash's second formal parameter 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'. * NEWS: Add NEWS entry. --- NEWS | 12 ++ doc/ref/api-data.texi | 44 ++++++-- libguile/numbers.c | 230 ++++++++++++++++++++++++++++------------- libguile/numbers.h | 3 +- test-suite/tests/numbers.test | 114 +++++++++------------ 5 files changed, 251 insertions(+), 152 deletions(-) diff --git a/NEWS b/NEWS index b53386a..c21e17c 100644 --- a/NEWS +++ b/NEWS @@ -1051,6 +1051,18 @@ R5RS integer-only operators `quotient' and `remainder'. For `round-quotient', `round-remainder', and `round/', Q is rounded to the nearest integer, with ties going to the nearest even integer. +*** New procedure: `round-ash', a rounding arithmetic shift operator + +round-ash is similar to ash, the arithmetic shift operator, except +that round-ash rounds to the nearest integer, with ties going to the +nearest even integer, whereas ash rounds toward negative infinity. + +Note that: (round-ash N COUNT) = round(N * 2^COUNT), +compared with: (ash N COUNT) = floor(N * 2^COUNT), + +except that round-ash and ash are computed more efficiently, and +require that N and COUNT be exact integers. + *** Complex number changes Guile is now able to represent non-real complex numbers whose diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi index 5bef926..8f7c35a 100644 --- a/doc/ref/api-data.texi +++ b/doc/ref/api-data.texi @@ -1659,19 +1659,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" @@ -1682,6 +1679,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 d0aacb7..4893790 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -4680,19 +4680,122 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0, } #undef FUNC_NAME +/* n must be an exact integer, and count > 0. + Returns n * 2^count. */ +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"); +} + +/* n must be an exact integer, and count > 0. + Returns floor(n / 2^count). */ +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"); +} + +/* n must be an exact integer, and count > 0. + Returns round(n / 2^count). */ +static SCM +round_right_shift_exact_integer (SCM n, long count) +{ + if (SCM_I_INUMP (n)) + { + scm_t_inum nn = SCM_I_INUM (n); + scm_t_inum qq = SCM_SRS (nn, count); + + if (count >= SCM_I_FIXNUM_BIT) + return SCM_INUM0; + else 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)) + { + if ((mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1) || + (mpz_odd_p (SCM_I_BIG_MPZ (q)))) + 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" - "\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" + (SCM n, SCM count), + "Return @math{floor(@var{n} * 2^@var{count})}, but computed\n" + "more efficiently. @var{n} and @var{count} must be exact\n" + "integers.\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" @@ -4703,79 +4806,58 @@ 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})}, but computed\n" + "more efficiently. @var{n} and @var{count} must be exact\n" + "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 ab96981..16d351e 100644 --- a/libguile/numbers.h +++ b/libguile/numbers.h @@ -204,7 +204,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 cb582ed..9a7d85b 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? ;;; @@ -4814,3 +4749,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.1 --=-=-=--