X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Fbignum.lisp;h=df1d328abc067851e437d80991f057dbaf84c5f8;hb=bfe145acc01eb7a43790173db4f08610ae9cb07a;hp=03f23377ab672d40899065e3c00ab4d75fc8a280;hpb=6cc71ab8ffad49f43895ad0a1df6885c81876687;p=sbcl.git

diff --git a/src/code/bignum.lisp b/src/code/bignum.lisp
index 03f2337..df1d328 100644
--- a/src/code/bignum.lisp
+++ b/src/code/bignum.lisp
@@ -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,7 +781,7 @@
 				       (%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