[PATCH] Improved log/log10, add round-ash

[Mark H Weaver <mhw@netris.org>]

[-- 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