;;; bignum-logical-and bignum-logical-ior bignum-logical-xor
;;; bignum-logical-not bignum-load-byte bignum-deposit-byte
;;; bignum-truncate bignum-plus-p bignum-compare make-small-bignum
-;;; bignum-logcount
+;;; bignum-logbitp bignum-logcount
;;; These symbols define the interface to the compiler:
;;; bignum-type bignum-element-type bignum-index %allocate-bignum
;;; %bignum-length %bignum-set-length %bignum-ref %bignum-set
\f
;;;; GCD
+;;; I'm not sure why I need this FTYPE declaration. Compiled by the
+;;; target compiler, it can deduce the return type fine, but without
+;;; it, we pay a heavy price in BIGNUM-GCD when compiled by the
+;;; cross-compiler. -- CSR, 2004-07-19
+(declaim (ftype (sfunction (bignum-type bignum-index bignum-type bignum-index)
+ (unsigned-byte 29))
+ bignum-factors-of-two))
+(defun bignum-factors-of-two (a len-a b len-b)
+ (declare (type bignum-index len-a len-b) (type bignum-type a b))
+ (do ((i 0 (1+ i))
+ (end (min len-a len-b)))
+ ((= i end) (error "Unexpected zero bignums?"))
+ (declare (type bignum-index i end))
+ (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
+ (unless (zerop or-digits)
+ (return (do ((j 0 (1+ j))
+ (or-digits or-digits (%ashr or-digits 1)))
+ ((oddp or-digits) (+ (* i digit-size) j))
+ (declare (type (mod 32) j))))))))
+
(defun bignum-gcd (a b)
- (declare (type bignum-type a b))
(let* ((a (if (%bignum-0-or-plusp a (%bignum-length a))
a
(negate-bignum a nil)))
(b (if (%bignum-0-or-plusp b (%bignum-length b))
b
- (negate-bignum b nil)))
- (len-a (%bignum-length a))
+ (negate-bignum b nil))))
+ (declare (type bignum-type a b))
+ (when (< a b)
+ (rotatef a b))
+ ;; While the length difference of A and B is sufficiently large,
+ ;; reduce using MOD (slowish, but it should equalize the sizes of
+ ;; A and B pretty quickly). After that, use the binary GCD
+ ;; algorithm to handle the rest. The initial reduction using MOD
+ ;; is sufficient to get rid of the embarrasing order of magnitude
+ ;; difference in GCD/LCM performance between SBCL and most other
+ ;; lisps.
+ ;;
+ ;; FIXME: Using a better algorithm (for example Weber's accelerated
+ ;; integer GCD) would be nice.
+ ;; -- JES, 2004-07-31
+ (loop until (and (= (%bignum-length b) 1) (zerop (%bignum-ref b 0))) do
+ (when (<= (%bignum-length a) (1+ (%bignum-length b)))
+ (return-from bignum-gcd (bignum-binary-gcd a b)))
+ (let ((rem (mod a b)))
+ (if (fixnump rem)
+ (setf a (make-small-bignum rem))
+ (setf a rem))
+ (rotatef a b)))
+ a))
+
+(defun bignum-binary-gcd (a b)
+ (declare (type bignum-type a b))
+ (let* ((len-a (%bignum-length a))
(len-b (%bignum-length b)))
- (declare (type bignum-index len-a len-b))
+ (declare (type bignum-index len-a len-b))
(with-bignum-buffers ((a-buffer len-a a)
(b-buffer len-b b)
(res-buffer (max len-a len-b)))
(bignum-buffer-ashift-right a len-a
(+ (* index digit-size)
increment)))))))
-
-(defun bignum-factors-of-two (a len-a b len-b)
- (declare (type bignum-index len-a len-b) (type bignum-type a))
- (do ((i 0 (1+ i))
- (end (min len-a len-b)))
- ((= i end) (error "Unexpected zero bignums?"))
- (declare (type bignum-index i end))
- (let ((or-digits (%logior (%bignum-ref a i) (%bignum-ref b i))))
- (unless (zerop or-digits)
- (return (do ((j 0 (1+ j))
- (or-digits or-digits (%ashr or-digits 1)))
- ((oddp or-digits) (+ (* i digit-size) j))
- (declare (type (mod 32) j))))))))
\f
;;;; negation
(%normalize-bignum res res-len))
res)))))
((> count bignum-len)
- 0)
+ (if (%bignum-0-or-plusp bignum bignum-len) 0 -1))
;; Since a FIXNUM should be big enough to address anything in
;; memory, including arrays of bits, and since arrays of bits
;; take up about the same space as corresponding fixnums, there
;; should be no way that we fall through to this case: any shift
;; right by a bignum should give zero. But let's check anyway:
- (t (error "bignum overflow: can't shift right by ~S")))))
+ (t (error "bignum overflow: can't shift right by ~S" count)))))
(defun bignum-ashift-right-digits (bignum digits)
(declare (type bignum-type bignum)
(t
(round-up))))))
\f
-;;;; integer length and logcount
+;;;; integer length and logbitp/logcount
(defun bignum-integer-length (bignum)
(declare (type bignum-type bignum))
(+ (integer-length (%fixnum-digit-with-correct-sign digit))
(* len-1 digit-size))))
+(defun bignum-logbitp (index bignum)
+ (declare (type bignum-type bignum))
+ (let ((len (%bignum-length bignum)))
+ (declare (type bignum-index len))
+ (multiple-value-bind (word-index bit-index)
+ (floor index digit-size)
+ (if (>= word-index len)
+ (not (bignum-plus-p bignum))
+ (not (zerop (logand (%bignum-ref bignum word-index)
+ (ash 1 bit-index))))))))
+
(defun bignum-logcount (bignum)
(declare (type bignum-type bignum))
(let* ((length (%bignum-length bignum))