* [PATCH] Improved log/log10, add round-ash
@ 2011-02-15 16:12 Mark H Weaver
2011-02-15 23:50 ` Ludovic Courtès
0 siblings, 1 reply; 3+ messages in thread
From: Mark H Weaver @ 2011-02-15 16:12 UTC (permalink / raw)
To: guile-devel
[-- Attachment #1: Type: text/plain, Size: 570 bytes --]
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
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Fix comment above number-theoretic division tests --]
[-- Type: text/x-diff, Size: 1059 bytes --]
From 4716bd257fb2d709c44d717ebb9ad2d8f4b50b8e Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
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
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: Improvements to `log' and `log10' --]
[-- Type: text/x-diff, Size: 10437 bytes --]
From 057df863b37bbc71010e647608a9daa7ae876afc Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
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+<PI>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
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #4: Add `round-ash', a rounding arithmetic shift operator --]
[-- Type: text/x-diff, Size: 18886 bytes --]
From 16ca1ed41ecb23a69fbd421f2f7d380bcdc82fd0 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
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
^ permalink raw reply related [flat|nested] 3+ messages in thread
* Re: [PATCH] Improved log/log10, add round-ash
2011-02-15 16:12 [PATCH] Improved log/log10, add round-ash Mark H Weaver
@ 2011-02-15 23:50 ` Ludovic Courtès
2011-02-16 0:39 ` Mark H Weaver
0 siblings, 1 reply; 3+ messages in thread
From: Ludovic Courtès @ 2011-02-15 23:50 UTC (permalink / raw)
To: guile-devel
Hi Mark,
Mark H Weaver <mhw@netris.org> writes:
> 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.
I applied them both. I’ll revert the latter tomorrow if it breaks on
some other platform.
The third one will be for later.
Thanks!
Ludo’.
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: [PATCH] Improved log/log10, add round-ash
2011-02-15 23:50 ` Ludovic Courtès
@ 2011-02-16 0:39 ` Mark H Weaver
0 siblings, 0 replies; 3+ messages in thread
From: Mark H Weaver @ 2011-02-16 0:39 UTC (permalink / raw)
To: Ludovic Courtès; +Cc: guile-devel
[-- Attachment #1: Type: text/plain, Size: 394 bytes --]
ludo@gnu.org (Ludovic Courtès) writes:
> I applied them both. I’ll revert the latter tomorrow if it breaks on
> some other platform.
I see one potential portability problem with the new patch.
Can you please apply the attached fix?
With this new patch applied, I think it should be portable,
assuming that the mysterious log1p gnulib module works.
Thanks,
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Portability fix for new log and log10 --]
[-- Type: text/x-diff, Size: 1027 bytes --]
From 873dfab119592658adce812a0fea0d6a688a3637 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Tue, 15 Feb 2011 19:29:41 -0500
Subject: [PATCH] Portability fix for new log and log10
* libguile/numbers.c: Define M_LN2 if it's not already defined.
Fix error in comment.
---
libguile/numbers.c | 5 ++++-
1 files changed, 4 insertions(+), 1 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index d0aacb7..b8cfa5d 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -72,6 +72,9 @@
#ifndef M_LOG10E
#define M_LOG10E 0.43429448190325182765
#endif
+#ifndef M_LN2
+#define M_LN2 0.69314718055994530942
+#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
@@ -9399,7 +9402,7 @@ log_of_exact_integer_with_size (SCM n, long size)
return log_of_shifted_double (scm_to_double (n), 0);
}
-/* Returns log(n), for exact integer n of integer-length size */
+/* Returns log(n), for exact integer n */
static SCM
log_of_exact_integer (SCM n)
{
--
1.7.1
^ permalink raw reply related [flat|nested] 3+ messages in thread
end of thread, other threads:[~2011-02-16 0:39 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2011-02-15 16:12 [PATCH] Improved log/log10, add round-ash Mark H Weaver
2011-02-15 23:50 ` Ludovic Courtès
2011-02-16 0:39 ` Mark H Weaver
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).