[PATCH] round-ash, a rounding arithmetic shift operator

[PATCH] round-ash, a rounding arithmetic shift operator

[Mark H Weaver]

[-- Attachment #1: Type: text/plain, Size: 557 bytes --]

The first patch is trivial, but there for the sake of correctness.
The second patch adds round-ash, a rounding arithmetic shift operator.

  (round-ash n count) ==> (round (* n (expt 2 count)))

but it's implemented much more efficiently than that, and requires that
both n and count are exact integers.  It cannot be implemented very
efficiently in scheme, and I needed it to normalize floating-point
significands using the default IEEE rounding mode.  I think it probably
has wider utility.  It would be great to have it in 2.0.  Any chance?

     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 24504319b81aeabd5ac7967ed4f69428b8a5fbea 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/2] 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: Add `round-ash', a rounding arithmetic shift operator --]
[-- Type: text/x-diff, Size: 18848 bytes --]

From d702946bcd1af2581e980387914b1086480564e6 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 14 Feb 2011 23:25:13 -0500
Subject: [PATCH 2/2] 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 6bebbf6..e473644 100644
--- a/NEWS
+++ b/NEWS
@@ -23,6 +23,18 @@ instead.
 `define-once' is like Lisp's `defvar': it creates a toplevel binding,
 but only if one does not exist already.
 
+** 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.
+
 ** Added four new sets of fast quotient and remainder operators
 
 Added four new sets of fast quotient and remainder operators with
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 59d8e74..a87b27b 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4682,19 +4682,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"
@@ -4705,79 +4808,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 9e9728f..7a75780 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?
 ;;;
 
@@ -4770,3 +4705,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

Re: [PATCH] round-ash, a rounding arithmetic shift operator

[Andy Wingo]

On Tue 15 Feb 2011 10:49, Mark H Weaver <mhw@netris.org> writes:

> The first patch is trivial, but there for the sake of correctness.

Please apply, thanks.

> The second patch adds round-ash, a rounding arithmetic shift operator.
>
>   (round-ash n count) ==> (round (* n (expt 2 count)))
>
> but it's implemented much more efficiently than that, and requires that
> both n and count are exact integers.  It cannot be implemented very
> efficiently in scheme, and I needed it to normalize floating-point
> significands using the default IEEE rounding mode.  I think it probably
> has wider utility.  It would be great to have it in 2.0.  Any chance?

I'm not sure I understand the name.  There is no need to call "round" on
the result of calling "ash".  Can you think of another name?

Cheers,

Andy
-- 
http://wingolog.org/

Re: [PATCH] round-ash, a rounding arithmetic shift operator

[Mark H Weaver]

Andy Wingo <wingo@pobox.com> writes:

> On Tue 15 Feb 2011 10:49, Mark H Weaver <mhw@netris.org> writes:
>> The first patch is trivial, but there for the sake of correctness.
>
> Please apply, thanks.

Ludo applied this before the 2.0.0 release.

>> The second patch adds round-ash, a rounding arithmetic shift operator.
>>
>>   (round-ash n count) ==> (round (* n (expt 2 count)))
>>
>> but it's implemented much more efficiently than that, and requires that
>> both n and count are exact integers.  It cannot be implemented very
>> efficiently in scheme, and I needed it to normalize floating-point
>> significands using the default IEEE rounding mode.  I think it probably
>> has wider utility.  It would be great to have it in 2.0.  Any chance?
>
> I'm not sure I understand the name.  There is no need to call "round" on
> the result of calling "ash".  Can you think of another name?

The name is inspired by Taylor Campbell's names for the division
operators.  There is no need to call "round" on the result of calling
"quotient", and yet we have adopted the name "round-quotient" for a
rounded quotient.  For consistency, I chose the name "round-ash" for a
rounded arithmetic shift.

For backward compatibility, we keep the name "quotient" which is short
for "truncate-quotient", and similarly "ash" is short for "floor-ash".

In theory, there could be an ash variant for every quotient operator,
where the divisor is constrained to be a power of 2.  However, the only
variants needed to implement the IEEE 754 rounding modes for the usual
sign-magnitude representation (where the significand is always
non-negative) are floor-ash, ceiling-ash, and round-ash:

  (define floor-ash ash)

  (define (ceiling-ash n count)
    (- (ash (- n) count)))

  (use-modules (srfi srfi-60))

  (define (round-ash n count)
    (let ((r (ash n count)))
      (if (and (negative? count)
               (bit-set? (- -1 count) n)
               (or (odd? r) (< (first-set-bit n) (- -1 count))))
          (1+ r)
          r)))

Given that round-ash is frequently called when using the default IEEE
754 rounding mode, it seems worth having a more efficient version than
the code above.

Having said all this, let's hold off on this patch for now.  I'm
modifying my bigfloat module to do strictly correct rounding, and that
may require something slightly more powerful than ash and round-ash.  In
particular, I it may be necessary to have variants that return the
remainder as well as the quotient, or maybe just some partial
information about the remainder.

     Thanks,
       Mark

Re: [PATCH] round-ash, a rounding arithmetic shift operator

[Andy Wingo]

On Tue 22 Feb 2011 18:54, Mark H Weaver <mhw@netris.org> writes:

> Andy Wingo <wingo@pobox.com> writes:
>
>> On Tue 15 Feb 2011 10:49, Mark H Weaver <mhw@netris.org> writes:
>>> The first patch is trivial, but there for the sake of correctness.
>>
>> Please apply, thanks.
>
> Ludo applied this before the 2.0.0 release.

Ah, cool.

>   (define (round-ash n count)
>     (let ((r (ash n count)))
>       (if (and (negative? count)
>                (bit-set? (- -1 count) n)
>                (or (odd? r) (< (first-set-bit n) (- -1 count))))
>           (1+ r)
>           r)))

Thanks for all the explanation.

You might get good results with (logand n (ash 1 X)), given that both
logand and ash have opcodes.  Dunno.

> Having said all this, let's hold off on this patch for now.  I'm
> modifying my bigfloat module to do strictly correct rounding, and that
> may require something slightly more powerful than ash and round-ash.  In
> particular, I it may be necessary to have variants that return the
> remainder as well as the quotient, or maybe just some partial
> information about the remainder.

OK, will hold off then.  Happy hacking!

Andy
-- 
http://wingolog.org/