summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMark H Weaver <mhw@netris.org>2013-03-03 04:35:09 -0500
committerMark H Weaver <mhw@netris.org>2013-03-12 15:39:20 -0400
commite08a12b5356c20ed0418bcaee136eb3632c5616f (patch)
treeac3a759997554ee2889ffa8e0d1a09f6ce3023ee
parenta285b18ca820e089e2e5d02f8ed07a1e341dffc3 (diff)
downloadguile-e08a12b5356c20ed0418bcaee136eb3632c5616f.tar.gz
Add 'round-ash', a rounding arithmetic shift operator
* libguile/numbers.c (left_shift_exact_integer, floor_right_shift_exact_integer, round_right_shift_exact_integer): New static functions. (scm_round_ash): New procedure. (scm_ash): Reimplement in terms of 'left_shift_exact_integer' and 'floor_right_shift_exact_integer'. * libguile/numbers.h: Add prototype for scm_round_ash. Rename the second argument of 'scm_ash' from 'cnt' to 'count'. * test-suite/tests/numbers.test (round-ash, ash): Add new unified testing framework for 'ash' and 'round-ash'. Previously, the tests for 'ash' were not very comprehensive; for example, they did not include a single test where the number to be shifted was a bignum. * doc/ref/api-data.texi (Bitwise Operations): Add documentation for 'round-ash'. Improve documentation for `ash'.
-rw-r--r--doc/ref/api-data.texi42
-rw-r--r--libguile/numbers.c226
-rw-r--r--libguile/numbers.h3
-rw-r--r--test-suite/tests/numbers.test114
4 files changed, 233 insertions, 152 deletions
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index fb12d2ccf..81c6d5b70 100644
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1686,19 +1686,15 @@ 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.
+@deffn {Scheme Procedure} ash n count
+@deffnx {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * 2^@var{count})}.
+@var{n} and @var{count} must be exact integers.
-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.)
-
-With @var{n} viewed as an infinite precision twos complement,
-@code{ash} means a left shift introducing zero bits, or a right shift
-dropping bits.
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{ash} means a left shift introducing zero bits
+when @var{count} is positive, or a right shift dropping bits
+when @var{count} is negative. This is an ``arithmetic'' shift.
@lisp
(number->string (ash #b1 3) 2) @result{} "1000"
@@ -1709,6 +1705,28 @@ 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})}.
+@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 2b64a748c..3f2afdebb 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4791,19 +4791,119 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
}
#undef FUNC_NAME
+/* Efficiently compute (N * 2^COUNT),
+ where N is an exact integer, and COUNT > 0. */
+static SCM
+left_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+
+ /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will always
+ overflow a non-zero fixnum. For smaller shifts we check the
+ bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
+ all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
+ Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - count)". */
+
+ if (nn == 0)
+ return n;
+ else if (count < SCM_I_FIXNUM_BIT-1 &&
+ ((scm_t_bits) (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - count)) + 1)
+ <= 1))
+ return SCM_I_MAKINUM (nn << count);
+ else
+ {
+ SCM result = scm_i_inum2big (nn);
+ mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
+ count);
+ return result;
+ }
+ }
+ else if (SCM_BIGP (n))
+ {
+ SCM result = scm_i_mkbig ();
+ mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n), count);
+ scm_remember_upto_here_1 (n);
+ return result;
+ }
+ else
+ scm_syserror ("left_shift_exact_integer");
+}
+
+/* Efficiently compute floor (N / 2^COUNT),
+ where N is an exact integer and COUNT > 0. */
+static SCM
+floor_right_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+
+ if (count >= SCM_I_FIXNUM_BIT)
+ return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM (-1));
+ else
+ return SCM_I_MAKINUM (SCM_SRS (nn, count));
+ }
+ else if (SCM_BIGP (n))
+ {
+ SCM result = scm_i_mkbig ();
+ mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
+ count);
+ scm_remember_upto_here_1 (n);
+ return scm_i_normbig (result);
+ }
+ else
+ scm_syserror ("floor_right_shift_exact_integer");
+}
+
+/* Efficiently compute round (N / 2^COUNT),
+ where N is an exact integer and COUNT > 0. */
+static SCM
+round_right_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ if (count >= SCM_I_FIXNUM_BIT)
+ return SCM_INUM0;
+ else
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+ scm_t_inum qq = SCM_SRS (nn, count);
+
+ if (0 == (nn & (1L << (count-1))))
+ return SCM_I_MAKINUM (qq); /* round down */
+ else if (nn & ((1L << (count-1)) - 1))
+ return SCM_I_MAKINUM (qq + 1); /* round up */
+ else
+ return SCM_I_MAKINUM ((~1L) & (qq + 1)); /* round to even */
+ }
+ }
+ else if (SCM_BIGP (n))
+ {
+ SCM q = scm_i_mkbig ();
+
+ mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), count);
+ if (mpz_tstbit (SCM_I_BIG_MPZ (n), count-1)
+ && (mpz_odd_p (SCM_I_BIG_MPZ (q))
+ || (mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1)))
+ mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
+ scm_remember_upto_here_1 (n);
+ return scm_i_normbig (q);
+ }
+ else
+ scm_syserror ("round_right_shift_exact_integer");
+}
+
SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
- (SCM n, SCM cnt),
- "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
- "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
+ (SCM n, SCM count),
+ "Return @math{floor(@var{n} * 2^@var{count})}.\n"
+ "@var{n} and @var{count} must be exact integers.\n"
"\n"
- "This is effectively a multiplication by 2^@var{cnt}, and when\n"
- "@var{cnt} is negative it's a division, rounded towards negative\n"
- "infinity. (Note that this is not the same rounding as\n"
- "@code{quotient} does.)\n"
- "\n"
- "With @var{n} viewed as an infinite precision twos complement,\n"
- "@code{ash} means a left shift introducing zero bits, or a right\n"
- "shift dropping bits.\n"
+ "With @var{n} viewed as an infinite-precision twos-complement\n"
+ "integer, @code{ash} means a left shift introducing zero bits\n"
+ "when @var{count} is positive, or a right shift dropping bits\n"
+ "when @var{count} is negative. This is an ``arithmetic'' shift.\n"
"\n"
"@lisp\n"
"(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
@@ -4814,79 +4914,57 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
"@end lisp")
#define FUNC_NAME s_scm_ash
{
- long bits_to_shift;
- bits_to_shift = scm_to_long (cnt);
-
- if (SCM_I_INUMP (n))
+ if (SCM_I_INUMP (n) || SCM_BIGP (n))
{
- scm_t_inum nn = SCM_I_INUM (n);
+ long bits_to_shift = scm_to_long (count);
if (bits_to_shift > 0)
- {
- /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
- overflow a non-zero fixnum. For smaller shifts we check the
- bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
- all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
- Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
- bits_to_shift)". */
-
- if (nn == 0)
- return n;
-
- if (bits_to_shift < SCM_I_FIXNUM_BIT-1
- && ((scm_t_bits)
- (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1)
- <= 1))
- {
- return SCM_I_MAKINUM (nn << bits_to_shift);
- }
- else
- {
- SCM result = scm_i_inum2big (nn);
- mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
- bits_to_shift);
- return result;
- }
- }
+ return left_shift_exact_integer (n, bits_to_shift);
+ else if (SCM_LIKELY (bits_to_shift < 0))
+ return floor_right_shift_exact_integer (n, -bits_to_shift);
else
- {
- bits_to_shift = -bits_to_shift;
- if (bits_to_shift >= SCM_LONG_BIT)
- return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM(-1));
- else
- return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift));
- }
-
+ return n;
}
- else if (SCM_BIGP (n))
- {
- SCM result;
+ else
+ SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
+}
+#undef FUNC_NAME
- if (bits_to_shift == 0)
- return n;
+SCM_DEFINE (scm_round_ash, "round-ash", 2, 0, 0,
+ (SCM n, SCM count),
+ "Return @math{round(@var{n} * 2^@var{count})}.\n"
+ "@var{n} and @var{count} must be exact integers.\n"
+ "\n"
+ "With @var{n} viewed as an infinite-precision twos-complement\n"
+ "integer, @code{round-ash} means a left shift introducing zero\n"
+ "bits when @var{count} is positive, or a right shift rounding\n"
+ "to the nearest integer (with ties going to the nearest even\n"
+ "integer) when @var{count} is negative. This is a rounded\n"
+ "``arithmetic'' shift.\n"
+ "\n"
+ "@lisp\n"
+ "(number->string (round-ash #b1 3) 2) @result{} \"1000\"\n"
+ "(number->string (round-ash #b1010 -1) 2) @result{} \"101\"\n"
+ "(number->string (round-ash #b1010 -2) 2) @result{} \"10\"\n"
+ "(number->string (round-ash #b1011 -2) 2) @result{} \"11\"\n"
+ "(number->string (round-ash #b1101 -2) 2) @result{} \"11\"\n"
+ "(number->string (round-ash #b1110 -2) 2) @result{} \"100\"\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_round_ash
+{
+ if (SCM_I_INUMP (n) || SCM_BIGP (n))
+ {
+ long bits_to_shift = scm_to_long (count);
- result = scm_i_mkbig ();
- if (bits_to_shift >= 0)
- {
- mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
- bits_to_shift);
- return result;
- }
+ if (bits_to_shift > 0)
+ return left_shift_exact_integer (n, bits_to_shift);
+ else if (SCM_LIKELY (bits_to_shift < 0))
+ return round_right_shift_exact_integer (n, -bits_to_shift);
else
- {
- /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
- we have to allocate a bignum even if the result is going to be a
- fixnum. */
- mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
- -bits_to_shift);
- return scm_i_normbig (result);
- }
-
+ return n;
}
else
- {
- SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
- }
+ SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
}
#undef FUNC_NAME
diff --git a/libguile/numbers.h b/libguile/numbers.h
index 2c8b2602e..912f287bb 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -206,7 +206,8 @@ SCM_API SCM scm_logbit_p (SCM n1, SCM n2);
SCM_API SCM scm_lognot (SCM n);
SCM_API SCM scm_modulo_expt (SCM n, SCM k, SCM m);
SCM_API SCM scm_integer_expt (SCM z1, SCM z2);
-SCM_API SCM scm_ash (SCM n, SCM cnt);
+SCM_API SCM scm_ash (SCM n, SCM count);
+SCM_API SCM scm_round_ash (SCM n, SCM count);
SCM_API SCM scm_bit_extract (SCM n, SCM start, SCM end);
SCM_API SCM scm_logcount (SCM n);
SCM_API SCM scm_integer_length (SCM n);
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index c4e819db2..bb1424853 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?
;;;
@@ -4914,3 +4849,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))