X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Fbignum.lisp;h=df1d328abc067851e437d80991f057dbaf84c5f8;hb=bfe145acc01eb7a43790173db4f08610ae9cb07a;hp=cf70586301619c9861adafa8ff315ad80bf5509c;hpb=93ba859423ec6e035a7b22a76a2ac70038691d65;p=sbcl.git diff --git a/src/code/bignum.lisp b/src/code/bignum.lisp index cf70586..df1d328 100644 --- a/src/code/bignum.lisp +++ b/src/code/bignum.lisp @@ -22,7 +22,7 @@ ;;; 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 @@ -519,17 +519,62 @@ ;;;; 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))) @@ -605,19 +650,6 @@ (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)))))))) ;;;; negation @@ -749,13 +781,13 @@ (%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) @@ -812,7 +844,7 @@ (bignum-ashift-left-unaligned bignum digits n-bits res-len)))) ;; Left shift by a number too big to be represented as a fixnum ;; would exceed our memory capacity, since a fixnum is big enough - ;; index any array, including a bit array. + ;; to index any array, including a bit array. (error "can't represent result of left shift"))) (defun bignum-ashift-left-digits (bignum bignum-len digits) @@ -1009,7 +1041,7 @@ (t (round-up)))))) -;;;; integer length and logcount +;;;; integer length and logbitp/logcount (defun bignum-integer-length (bignum) (declare (type bignum-type bignum)) @@ -1021,6 +1053,17 @@ (+ (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))