X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Fbignum.lisp;h=d05a6fadf6414d9069f6a330e584d5f98b61de8d;hb=77d94d36bcfd3d5eea73ad51e6ee621a8938f995;hp=e50ec3270701c57ec74c0707fae440f6e6c2131c;hpb=b5dc433da5b8bd3b36db88f7ec35cdb03cb64384;p=sbcl.git

diff --git a/src/code/bignum.lisp b/src/code/bignum.lisp
index e50ec32..d05a6fa 100644
--- a/src/code/bignum.lisp
+++ b/src/code/bignum.lisp
@@ -112,6 +112,8 @@
 
 (defconstant maximum-bignum-length (1- (ash 1 (- sb!vm:n-word-bits
                                                  sb!vm:n-widetag-bits))))
+
+(defconstant all-ones-digit (1- (ash 1 sb!vm:n-word-bits)))
 
 ;;;; internal inline routines
 
@@ -374,12 +376,17 @@
 ;;; Operations requiring a subtraction without the overhead of intermediate
 ;;; results, such as GCD, use this. It assumes Result is big enough for the
 ;;; result.
-(defun subtract-bignum-buffers (a len-a b len-b result)
+(defun subtract-bignum-buffers-with-len (a len-a b len-b result len-res)
   (declare (type bignum-type a b)
 	   (type bignum-index len-a len-b))
-  (let ((len-res (max len-a len-b)))
-    (subtract-bignum-loop a len-a b len-b result len-res
-			  %normalize-bignum-buffer)))
+  (subtract-bignum-loop a len-a b len-b result len-res
+			%normalize-bignum-buffer))
+
+(defun subtract-bignum-buffers (a len-a b len-b result)
+  (declare (type bignum-type a b)
+	   (type bignum-index len-a len-b))  
+  (subtract-bignum-loop a len-a b len-b result (max len-a len-b)
+			%normalize-bignum-buffer))
 
 ;;;; multiplication
 
@@ -521,6 +528,16 @@
 
 ;;;; GCD
 
+(eval-when (:compile-toplevel :load-toplevel :execute)
+  ;; The asserts in the GCD implementation are way too expensive to
+  ;; check in normal use, and are disabled here.
+  (sb!xc:defmacro gcd-assert (&rest args)
+    (if nil
+	`(assert ,@args)))
+  ;; We'll be doing a lot of modular arithmetic.
+  (sb!xc:defmacro M (form)
+    `(logand all-ones-digit ,form)))
+
 ;;; 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
@@ -541,6 +558,127 @@
 		    ((oddp or-digits) (+ (* i digit-size) j))
 		  (declare (type (mod #.sb!vm:n-word-bits) j))))))))
 
+;;; Multiply a bignum buffer with a fixnum or a digit, storing the
+;;; result in another bignum buffer, and without using any
+;;; temporaries. Inlined to avoid boxing smallnum if it's actually a
+;;; digit. Needed by GCD, should possibly OAOO with
+;;; MULTIPLY-BIGNUM-AND-FIXNUM.
+(declaim (inline multiply-bignum-buffer-and-smallnum-to-buffer))
+(defun multiply-bignum-buffer-and-smallnum-to-buffer (bignum bignum-len
+							     smallnum res)
+  (declare (type bignum-type bignum))
+  (let* ((bignum-plus-p (%bignum-0-or-plusp bignum bignum-len))
+	 (smallnum-plus-p (not (minusp smallnum)))
+	 (smallnum (if smallnum-plus-p smallnum (- smallnum)))
+	 (carry-digit 0)) 
+    (declare (type bignum-type bignum res)
+	     (type bignum-index bignum-len)
+	     (type bignum-element-type smallnum carry-digit))
+    (unless bignum-plus-p
+      (negate-bignum-buffer-in-place bignum bignum-len))
+    (dotimes (index bignum-len)
+      (declare (type bignum-index index))
+      (multiple-value-bind (next-digit low)
+	  (%multiply-and-add (%bignum-ref bignum index)
+			     smallnum
+			     carry-digit)
+	(declare (type bignum-element-type next-digit low))
+	(setf carry-digit next-digit)
+	(setf (%bignum-ref res index) low)))
+    (setf (%bignum-ref res bignum-len) carry-digit)
+    (unless bignum-plus-p
+      (negate-bignum-buffer-in-place bignum bignum-len))
+    (let ((res-len (%normalize-bignum-buffer res (1+ bignum-len))))
+      (unless (eq bignum-plus-p smallnum-plus-p)
+	(negate-bignum-buffer-in-place res res-len))
+      res-len)))
+
+;;; Given U and V, return U / V mod 2^32. Implements the algorithm in the
+;;; paper, but uses some clever bit-twiddling nicked from Nickle to do it.
+(declaim (inline bmod))
+(defun bmod (u v)
+  (let ((ud (%bignum-ref u 0))
+	(vd (%bignum-ref v 0))
+	(umask 0)
+	(imask 1)
+	(m 0))
+    (declare (type (unsigned-byte #.sb!vm:n-word-bits) ud vd umask imask m))
+    (dotimes (i digit-size)
+      (setf umask (logior umask imask))
+      (unless (zerop (logand ud umask))
+	(setf ud (M (- ud vd)))
+	(setf m (M (logior m imask))))
+      (setf imask (M (ash imask 1)))
+      (setf vd (M (ash vd 1))))
+    m))
+
+(defun dmod (u u-len v v-len tmp1)
+  (loop while (> (bignum-buffer-integer-length u u-len)
+		 (+ (bignum-buffer-integer-length v v-len)
+		    digit-size))
+    do
+    (unless (zerop (%bignum-ref u 0))
+      (let* ((bmod (bmod u v))
+	     (tmp1-len (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
+								      bmod
+								      tmp1)))
+	(setf u-len (subtract-bignum-buffers u u-len
+					     tmp1 tmp1-len
+					     u))
+	(bignum-abs-buffer u u-len)))
+    (gcd-assert (zerop (%bignum-ref u 0)))
+    (setf u-len (bignum-buffer-ashift-right u u-len digit-size)))
+  (let* ((d (+ 1 (- (bignum-buffer-integer-length u u-len)
+		    (bignum-buffer-integer-length v v-len))))
+	 (n (1- (ash 1 d))))
+    (declare (type (unsigned-byte #.(integer-length #.sb!vm:n-word-bits)) d)
+	     (type (unsigned-byte #.sb!vm:n-word-bits) n))
+    (gcd-assert (>= d 0))
+    (unless (zerop (logand (%bignum-ref u 0) n))
+      (let ((tmp1-len
+	     (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
+							    (logand n (bmod u
+									    v))
+							    tmp1)))
+	(setf u-len (subtract-bignum-buffers u u-len
+					     tmp1 tmp1-len
+					     u))
+	(bignum-abs-buffer u u-len)))
+    u-len))
+
+(defconstant lower-ones-digit (1- (ash 1 (truncate sb!vm:n-word-bits 2))))
+  
+;;; Find D and N such that (LOGAND ALL-ONES-DIGIT (- (* D X) (* N Y))) is 0,
+;;; (< 0 N LOWER-ONES-DIGIT) and (< 0 (ABS D) LOWER-ONES-DIGIT).
+(defun reduced-ratio-mod (x y)
+  (let* ((c (bmod x y))
+	 (n1 c)
+	 (d1 1)
+	 (n2 (M (1+ (M (lognot n1)))))
+	 (d2 (M -1)))
+    (declare (type (unsigned-byte #.sb!vm:n-word-bits) n1 d1 n2 d2))
+    (loop while (> n2 (expt 2 (truncate digit-size 2))) do
+	  (loop for i of-type (mod #.sb!vm:n-word-bits)
+		downfrom (- (integer-length n1) (integer-length n2))
+		while (>= n1 n2) do
+		(when (>= n1 (M (ash n2 i)))
+		  (psetf n1 (M (- n1 (M (ash n2 i))))
+			 d1 (M (- d1 (M (ash d2 i)))))))
+	  (psetf n1 n2
+		 d1 d2
+		 n2 n1
+		 d2 d1))    
+    (values n2 (if (>= d2 (expt 2 (1- digit-size)))
+		   (lognot (logand most-positive-fixnum (lognot d2)))
+		   (logand lower-ones-digit d2)))))
+
+
+(defun copy-bignum (a &optional (len (%bignum-length a)))
+  (let ((b (%allocate-bignum len)))
+    (bignum-replace b a)
+    (%bignum-set-length b len)
+    b))
+    
 ;;; Allocate a single word bignum that holds fixnum. This is useful when
 ;;; we are trying to mix fixnum and bignum operands.
 #!-sb-fluid (declaim (inline make-small-bignum))
@@ -549,39 +687,131 @@
     (setf (%bignum-ref res 0) (%fixnum-to-digit fixnum))
     res))
 
-(defun bignum-gcd (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))))
-    (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)))
-    (if (= (%bignum-length a) 1)
-	(%normalize-bignum a 1)
-	a)))
-  
+;; When the larger number is less than this many bignum digits long, revert
+;; to old algorithm. 
+(defparameter *accelerated-gcd-cutoff* 3)
+
+;;; Alternate between k-ary reduction with the help of
+;;; REDUCED-RATIO-MOD and digit modulus reduction via DMOD. Once the
+;;; arguments get small enough, drop through to BIGNUM-MOD-GCD (since
+;;; k-ary reduction can introduce spurious factors, which need to be
+;;; filtered out). Reference: Kenneth Weber, "The accelerated integer
+;;; GCD algorithm", ACM Transactions on Mathematical Software, volume
+;;; 21, number 1, March 1995, epp. 111-122.
+(defun bignum-gcd (u0 v0)
+  (declare (type bignum-type u0 v0))
+  (let* ((u1 (if (%bignum-0-or-plusp u0 (%bignum-length u0))
+		 u0
+		 (negate-bignum u0 nil)))
+	 (v1 (if (%bignum-0-or-plusp v0 (%bignum-length v0))
+		 v0
+		 (negate-bignum v0 nil))))
+    (if (zerop v1)
+	(return-from bignum-gcd u1))	
+    (when (> u1 v1)
+      (rotatef u1 v1))
+    (let ((n (mod v1 u1)))
+      (setf v1 (if (fixnump n)
+		   (make-small-bignum n)
+		   n)))
+    (if (and (= 1 (%bignum-length v1))
+	     (zerop (%bignum-ref v1 0)))
+	(return-from bignum-gcd (%normalize-bignum u1
+						   (%bignum-length u1))))
+    (let* ((buffer-len (+ 2 (%bignum-length u1)))
+	   (u (%allocate-bignum buffer-len))
+	   (u-len (%bignum-length u1))
+	   (v (%allocate-bignum buffer-len))
+	   (v-len (%bignum-length v1))
+	   (tmp1 (%allocate-bignum buffer-len))
+	   (tmp1-len 0)
+	   (tmp2 (%allocate-bignum buffer-len))
+	   (tmp2-len 0)
+	   (factors-of-two
+	    (bignum-factors-of-two u1 (%bignum-length u1)
+				   v1 (%bignum-length v1))))
+      (declare (type (or null bignum-index)
+		     buffer-len u-len v-len tmp1-len tmp2-len))
+      (bignum-replace u u1)
+      (bignum-replace v v1)
+      (setf u-len
+	    (make-gcd-bignum-odd u
+				 (bignum-buffer-ashift-right u u-len
+							     factors-of-two)))
+      (setf v-len
+	    (make-gcd-bignum-odd v
+				 (bignum-buffer-ashift-right v v-len
+							     factors-of-two)))
+      (loop until (or (< u-len *accelerated-gcd-cutoff*)
+		      (not v-len)
+		      (zerop v-len)
+		      (and (= 1 v-len)
+			   (zerop (%bignum-ref v 0))))	
+	do
+	(gcd-assert (= buffer-len (%bignum-length u)
+		       (%bignum-length v)
+		       (%bignum-length tmp1)
+		       (%bignum-length tmp2)))
+	(if (> (bignum-buffer-integer-length u u-len)
+	       (+ #.(truncate sb!vm:n-word-bits 4)
+		  (bignum-buffer-integer-length v v-len)))
+	    (setf u-len (dmod u u-len
+			      v v-len
+			      tmp1))
+	    (multiple-value-bind (n d) (reduced-ratio-mod u v)
+	      (setf tmp1-len
+		    (multiply-bignum-buffer-and-smallnum-to-buffer v v-len
+								   n tmp1))
+	      (setf tmp2-len
+		    (multiply-bignum-buffer-and-smallnum-to-buffer u u-len
+								   d tmp2))
+	      (gcd-assert (= (copy-bignum tmp2 tmp2-len)
+			     (* (copy-bignum u u-len) d)))
+	      (gcd-assert (= (copy-bignum tmp1 tmp1-len)
+			     (* (copy-bignum v v-len) n)))
+	      (setf u-len
+		    (subtract-bignum-buffers-with-len tmp1 tmp1-len
+						      tmp2 tmp2-len
+						      u
+						      (1+ (max tmp1-len
+							       tmp2-len))))
+	      (gcd-assert (or (zerop (- (copy-bignum tmp1 tmp1-len)
+					(copy-bignum tmp2 tmp2-len)))
+			      (= (copy-bignum u u-len)
+				 (- (copy-bignum tmp1 tmp1-len)
+				    (copy-bignum tmp2 tmp2-len)))))
+	      (bignum-abs-buffer u u-len)
+	      (gcd-assert (zerop (M u)))))
+	(setf u-len (make-gcd-bignum-odd u u-len))
+	(rotatef u v)	
+	(rotatef u-len v-len))
+      (setf u (copy-bignum u u-len))
+      (let ((n (bignum-mod-gcd v1 u)))	
+	(ash (bignum-mod-gcd u1 (if (fixnump n)
+				    (make-small-bignum n)
+				    n))
+	     factors-of-two)))))
+
+(defun bignum-mod-gcd (a b)
+  (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.
+  (loop until (and (= (%bignum-length b) 1) (zerop (%bignum-ref b 0))) do
+	(when (<= (%bignum-length a) (1+ (%bignum-length b)))
+	  (return-from bignum-mod-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)))
+  (if (= (%bignum-length a) 1)
+      (%normalize-bignum a 1)
+      a))
+
 (defun bignum-binary-gcd (a b)
   (declare (type bignum-type a b))
   (let* ((len-a (%bignum-length a))
@@ -635,7 +865,8 @@
 		   ((%digit-greater a-digit b-digit)
 		    (return
 		     (values a b len-b res
-			     (subtract-bignum-buffers a len-a b len-b res))))
+			     (subtract-bignum-buffers a len-a b len-b
+						      res))))
 		   (t
 		    (return
 		     (values b a len-a res
@@ -662,6 +893,7 @@
 		     (bignum-buffer-ashift-right a len-a
 						 (+ (* index digit-size)
 						    increment)))))))
+
 
 ;;;; negation
 
@@ -718,14 +950,20 @@
 
 ;;; This assumes bignum is positive; that is, the result of negating it will
 ;;; stay in the provided allocated bignum.
-(defun negate-bignum-in-place (bignum)
-  (bignum-negate-loop bignum (%bignum-length bignum) bignum)
+(defun negate-bignum-buffer-in-place (bignum bignum-len)
+  (bignum-negate-loop bignum bignum-len bignum)
   bignum)
+
+(defun negate-bignum-in-place (bignum)
+  (declare (inline negate-bignum-buffer-in-place))
+  (negate-bignum-buffer-in-place bignum (%bignum-length bignum)))
+
+(defun bignum-abs-buffer (bignum len)
+  (unless (%bignum-0-or-plusp bignum len)
+    (negate-bignum-in-place bignum len)))
 
 ;;;; shifting
 
-(defconstant all-ones-digit (1- (ash 1 sb!vm:n-word-bits)))
-
 (eval-when (:compile-toplevel :execute)
 
 ;;; This macro is used by BIGNUM-ASHIFT-RIGHT, BIGNUM-BUFFER-ASHIFT-RIGHT, and
@@ -1055,16 +1293,19 @@
 
 ;;;; integer length and logbitp/logcount
 
-(defun bignum-integer-length (bignum)
+(defun bignum-buffer-integer-length (bignum len)
   (declare (type bignum-type bignum))
-  (let* ((len (%bignum-length bignum))
-	 (len-1 (1- len))
+  (let* ((len-1 (1- len))
 	 (digit (%bignum-ref bignum len-1)))
     (declare (type bignum-index len len-1)
 	     (type bignum-element-type digit))
     (+ (integer-length (%fixnum-digit-with-correct-sign digit))
        (* len-1 digit-size))))
 
+(defun bignum-integer-length (bignum)
+  (declare (type bignum-type bignum))
+  (bignum-buffer-integer-length bignum (%bignum-length bignum)))
+
 (defun bignum-logbitp (index bignum)
   (declare (type bignum-type bignum))
   (let ((len (%bignum-length bignum)))
@@ -1761,274 +2002,303 @@ IS LESS EFFICIENT BUT EASIER TO MAINTAIN. BILL SAYS THIS CODE CERTAINLY WORKS!
 ;;;
 ;;; Normalize quotient and remainder. Cons result if necessary.
 
-;;; These are used by BIGNUM-TRUNCATE and friends in the general case.
-(defvar *truncate-x*)
-(defvar *truncate-y*)
-
-;;; Divide X by Y returning the quotient and remainder. In the
-;;; general case, we shift Y to set up for the algorithm, and we use
-;;; two buffers to save consing intermediate values. X gets
-;;; destructively modified to become the remainder, and we have to
-;;; shift it to account for the initial Y shift. After we multiple
-;;; bind q and r, we first fix up the signs and then return the
-;;; normalized results.
+
+;;; This used to be split into multiple functions, which shared state
+;;; in special variables *TRUNCATE-X* and *TRUNCATE-Y*. Having so many
+;;; special variable accesses in tight inner loops was having a large
+;;; effect on performance, so the helper functions have now been
+;;; refactored into local functions and the special variables into
+;;; lexicals.  There was also a lot of boxing and unboxing of
+;;; (UNSIGNED-BYTE 32)'s going on, which this refactoring
+;;; eliminated. This improves the performance on some CL-BENCH tests
+;;; by up to 50%, which is probably signigicant enough to justify the
+;;; reduction in readability that was introduced. --JES, 2004-08-07
 (defun bignum-truncate (x y)
   (declare (type bignum-type x y))
-  (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
-	 (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
-	 (x (if x-plusp x (negate-bignum x nil)))
-	 (y (if y-plusp y (negate-bignum y nil)))
-	 (len-x (%bignum-length x))
-	 (len-y (%bignum-length y)))
-    (multiple-value-bind (q r)
-	(cond ((< len-y 2)
-	       (bignum-truncate-single-digit x len-x y))
-	      ((plusp (bignum-compare y x))
-	       (let ((res (%allocate-bignum len-x)))
-		 (dotimes (i len-x)
-		   (setf (%bignum-ref res i) (%bignum-ref x i)))
-		 (values 0 res)))
-	      (t
-	       (let ((len-x+1 (1+ len-x)))
-		 (with-bignum-buffers ((*truncate-x* len-x+1)
-				       (*truncate-y* (1+ len-y)))
-		   (let ((y-shift (shift-y-for-truncate y)))
-		     (shift-and-store-truncate-buffers x len-x y len-y y-shift)
-		     (values (return-quotient-leaving-remainder len-x+1 len-y)
-			     ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
-			     ;; has executed, we just tidy up the remainder
-			     ;; (in *TRUNCATE-X*) and return it.
-			     (cond
-			      ((zerop y-shift)
-			       (let ((res (%allocate-bignum len-y)))
-				 (declare (type bignum-type res))
-				 (bignum-replace res *truncate-x* :end2 len-y)
-				 (%normalize-bignum res len-y)))
-			      (t
-			       (shift-right-unaligned
-				*truncate-x* 0 y-shift len-y
-				((= j res-len-1)
-				 (setf (%bignum-ref res j)
-				       (%ashr (%bignum-ref *truncate-x* i)
-					      y-shift))
-				 (%normalize-bignum res res-len))
-				res)))))))))
-      (let ((quotient (cond ((eq x-plusp y-plusp) q)
-			    ((typep q 'fixnum) (the fixnum (- q)))
-			    (t (negate-bignum-in-place q))))
-	    (rem (cond (x-plusp r)
-		       ((typep r 'fixnum) (the fixnum (- r)))
-		       (t (negate-bignum-in-place r)))))
-	(values (if (typep quotient 'fixnum)
-		    quotient
-		    (%normalize-bignum quotient (%bignum-length quotient)))
-		(if (typep rem 'fixnum)
-		    rem
-		    (%normalize-bignum rem (%bignum-length rem))))))))
-
-;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
-;;; fixes up the quotient and remainder with respect to sign and
-;;; normalization.
-;;;
-;;; We don't have to worry about shifting Y to make its most
-;;; significant digit sufficiently large for %FLOOR to return digit-size
-;;; quantities for the q-digit and r-digit. If Y is a single digit
-;;; bignum, it is already large enough for %FLOOR. That is, it has
-;;; some bits on pretty high in the digit.
-(defun bignum-truncate-single-digit (x len-x y)
-  (declare (type bignum-index len-x))
-  (let ((q (%allocate-bignum len-x))
-	(r 0)
-	(y (%bignum-ref y 0)))
-    (declare (type bignum-element-type r y))
-    (do ((i (1- len-x) (1- i)))
-	((minusp i))
-      (multiple-value-bind (q-digit r-digit) (%floor r (%bignum-ref x i) y)
-	(declare (type bignum-element-type q-digit r-digit))
-	(setf (%bignum-ref q i) q-digit)
-	(setf r r-digit)))
-    (let ((rem (%allocate-bignum 1)))
-      (setf (%bignum-ref rem 0) r)
-      (values q rem))))
-
-;;; a helper function for BIGNUM-TRUNCATE
-;;;
-;;; Divide *TRUNCATE-X* by *TRUNCATE-Y*, returning the quotient
-;;; and destructively modifying *TRUNCATE-X* so that it holds
-;;; the remainder.
-;;;
-;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
-;;;
-;;; *TRUNCATE-X* definitely has at least three digits, and it has one
-;;; more than *TRUNCATE-Y*. This keeps i, i-1, i-2, and low-x-digit
-;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
-(defun return-quotient-leaving-remainder (len-x len-y)
-  (declare (type bignum-index len-x len-y))
-  (let* ((len-q (- len-x len-y))
-	 ;; Add one for extra sign digit in case high bit is on.
-	 (q (%allocate-bignum (1+ len-q)))
-	 (k (1- len-q))
-	 (y1 (%bignum-ref *truncate-y* (1- len-y)))
-	 (y2 (%bignum-ref *truncate-y* (- len-y 2)))
-	 (i (1- len-x))
-	 (i-1 (1- i))
-	 (i-2 (1- i-1))
-	 (low-x-digit (- i len-y)))
-    (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
-	     (type bignum-element-type y1 y2))
-    (loop
-      (setf (%bignum-ref q k)
-	    (try-bignum-truncate-guess
-	     ;; This modifies *TRUNCATE-X*. Must access elements each pass.
-	     (bignum-truncate-guess y1 y2
-				    (%bignum-ref *truncate-x* i)
-				    (%bignum-ref *truncate-x* i-1)
-				    (%bignum-ref *truncate-x* i-2))
-	     len-y low-x-digit))
-      (cond ((zerop k) (return))
-	    (t (decf k)
-	       (decf low-x-digit)
-	       (shiftf i i-1 i-2 (1- i-2)))))
-    q))
-
-;;; This takes a digit guess, multiplies it by *TRUNCATE-Y* for a
-;;; result one greater in length than LEN-Y, and subtracts this result
-;;; from *TRUNCATE-X*. LOW-X-DIGIT is the first digit of X to start
-;;; the subtraction, and we know X is long enough to subtract a LEN-Y
-;;; plus one length bignum from it. Next we check the result of the
-;;; subtraction, and if the high digit in X became negative, then our
-;;; guess was one too big. In this case, return one less than GUESS
-;;; passed in, and add one value of Y back into X to account for
-;;; subtracting one too many. Knuth shows that the guess is wrong on
-;;; the order of 3/b, where b is the base (2 to the digit-size power)
-;;; -- pretty rarely.
-(defun try-bignum-truncate-guess (guess len-y low-x-digit)
-  (declare (type bignum-index low-x-digit len-y)
-	   (type bignum-element-type guess))
-  (let ((carry-digit 0)
-	(borrow 1)
-	(i low-x-digit))
-    (declare (type bignum-element-type carry-digit)
-	     (type bignum-index i)
-	     (fixnum borrow))
-    ;; Multiply guess and divisor, subtracting from dividend simultaneously.
-    (dotimes (j len-y)
-      (multiple-value-bind (high-digit low-digit)
-	  (%multiply-and-add guess
-			     (%bignum-ref *truncate-y* j)
-			     carry-digit)
-	(declare (type bignum-element-type high-digit low-digit))
-	(setf carry-digit high-digit)
-	(multiple-value-bind (x temp-borrow)
-	    (%subtract-with-borrow (%bignum-ref *truncate-x* i)
-				   low-digit
-				   borrow)
-	  (declare (type bignum-element-type x)
-		   (fixnum temp-borrow))
-	  (setf (%bignum-ref *truncate-x* i) x)
-	  (setf borrow temp-borrow)))
-      (incf i))
-    (setf (%bignum-ref *truncate-x* i)
-	  (%subtract-with-borrow (%bignum-ref *truncate-x* i)
-				 carry-digit borrow))
-    ;; See whether guess is off by one, adding one Y back in if necessary.
-    (cond ((%digit-0-or-plusp (%bignum-ref *truncate-x* i))
-	   guess)
-	  (t
-	   ;; If subtraction has negative result, add one divisor value back
-	   ;; in. The guess was one too large in magnitude.
-	   (let ((i low-x-digit)
-		 (carry 0))
+  (let (truncate-x truncate-y)
+    (labels	       
+        ;;; Divide X by Y when Y is a single bignum digit. BIGNUM-TRUNCATE
+        ;;; fixes up the quotient and remainder with respect to sign and
+        ;;; normalization.
+	;;;
+	;;; We don't have to worry about shifting Y to make its most
+	;;; significant digit sufficiently large for %FLOOR to return
+	;;; digit-size quantities for the q-digit and r-digit. If Y is
+	;;; a single digit bignum, it is already large enough for
+	;;; %FLOOR. That is, it has some bits on pretty high in the
+	;;; digit.
+	((bignum-truncate-single-digit (x len-x y)
+	   (declare (type bignum-index len-x))
+	   (let ((q (%allocate-bignum len-x))
+		 (r 0)
+		 (y (%bignum-ref y 0)))
+	     (declare (type bignum-element-type r y))
+	     (do ((i (1- len-x) (1- i)))
+		 ((minusp i))
+	       (multiple-value-bind (q-digit r-digit)
+		   (%floor r (%bignum-ref x i) y)
+		 (declare (type bignum-element-type q-digit r-digit))
+		 (setf (%bignum-ref q i) q-digit)
+		 (setf r r-digit)))
+	     (let ((rem (%allocate-bignum 1)))
+	       (setf (%bignum-ref rem 0) r)
+	       (values q rem))))	 
+	;;; This returns a guess for the next division step. Y1 is the
+	;;; highest y digit, and y2 is the second to highest y
+	;;; digit. The x... variables are the three highest x digits
+	;;; for the next division step.
+	;;;
+	;;; From Knuth, our guess is either all ones or x-i and x-i-1
+	;;; divided by y1, depending on whether x-i and y1 are the
+	;;; same. We test this guess by determining whether guess*y2
+	;;; is greater than the three high digits of x minus guess*y1
+	;;; shifted left one digit:
+	;;;    ------------------------------
+	;;;   |    x-i    |   x-i-1  | x-i-2 |
+	;;;    ------------------------------
+	;;;    ------------------------------
+	;;; - | g*y1 high | g*y1 low |   0   |
+	;;;    ------------------------------
+	;;;		...		  <   guess*y2     ???	 
+	;;; If guess*y2 is greater, then we decrement our guess by one
+	;;; and try again.  This returns a guess that is either
+	;;; correct or one too large.
+	 (bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
+	   (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
+	   (let ((guess (if (%digit-compare x-i y1)
+			    all-ones-digit
+			    (%floor x-i x-i-1 y1))))
+	     (declare (type bignum-element-type guess))
+	     (loop
+		 (multiple-value-bind (high-guess*y1 low-guess*y1)
+		     (%multiply guess y1)
+		   (declare (type bignum-element-type low-guess*y1
+				  high-guess*y1))
+		   (multiple-value-bind (high-guess*y2 low-guess*y2)
+		       (%multiply guess y2)
+		     (declare (type bignum-element-type high-guess*y2
+				    low-guess*y2))
+		     (multiple-value-bind (middle-digit borrow)
+			 (%subtract-with-borrow x-i-1 low-guess*y1 1)
+		       (declare (type bignum-element-type middle-digit)
+				(fixnum borrow))
+		       ;; Supplying borrow of 1 means there was no
+		       ;; borrow, and we know x-i-2 minus 0 requires
+		       ;; no borrow.
+		       (let ((high-digit (%subtract-with-borrow x-i
+								high-guess*y1
+								borrow)))
+			 (declare (type bignum-element-type high-digit))
+			 (if (and (%digit-compare high-digit 0)
+				  (or (%digit-greater high-guess*y2
+						      middle-digit)
+				      (and (%digit-compare middle-digit
+							   high-guess*y2)
+					   (%digit-greater low-guess*y2
+							   x-i-2))))
+			     (setf guess (%subtract-with-borrow guess 1 1))
+			     (return guess)))))))))
+	;;; Divide TRUNCATE-X by TRUNCATE-Y, returning the quotient
+	;;; and destructively modifying TRUNCATE-X so that it holds
+	;;; the remainder.
+	;;;
+	;;; LEN-X and LEN-Y tell us how much of the buffers we care about.
+	;;;
+	;;; TRUNCATE-X definitely has at least three digits, and it has one
+	;;; more than TRUNCATE-Y. This keeps i, i-1, i-2, and low-x-digit
+	;;; happy. Thanks to SHIFT-AND-STORE-TRUNCATE-BUFFERS.
+	 (return-quotient-leaving-remainder (len-x len-y)
+	   (declare (type bignum-index len-x len-y))
+	   (let* ((len-q (- len-x len-y))
+		  ;; Add one for extra sign digit in case high bit is on.
+		  (q (%allocate-bignum (1+ len-q)))
+		  (k (1- len-q))
+		  (y1 (%bignum-ref truncate-y (1- len-y)))
+		  (y2 (%bignum-ref truncate-y (- len-y 2)))
+		  (i (1- len-x))
+		  (i-1 (1- i))
+		  (i-2 (1- i-1))
+		  (low-x-digit (- i len-y)))
+	     (declare (type bignum-index len-q k i i-1 i-2 low-x-digit)
+		      (type bignum-element-type y1 y2))
+	     (loop
+		 (setf (%bignum-ref q k)
+		       (try-bignum-truncate-guess
+			;; This modifies TRUNCATE-X. Must access
+			;; elements each pass.
+			(bignum-truncate-guess y1 y2
+					       (%bignum-ref truncate-x i)
+					       (%bignum-ref truncate-x i-1)
+					       (%bignum-ref truncate-x i-2))
+			len-y low-x-digit))
+		 (cond ((zerop k) (return))
+		       (t (decf k)
+			  (decf low-x-digit)
+			  (shiftf i i-1 i-2 (1- i-2)))))
+	     q))
+	;;; This takes a digit guess, multiplies it by TRUNCATE-Y for a
+	;;; result one greater in length than LEN-Y, and subtracts this result
+	;;; from TRUNCATE-X. LOW-X-DIGIT is the first digit of X to start
+	;;; the subtraction, and we know X is long enough to subtract a LEN-Y
+	;;; plus one length bignum from it. Next we check the result of the
+	;;; subtraction, and if the high digit in X became negative, then our
+	;;; guess was one too big. In this case, return one less than GUESS
+	;;; passed in, and add one value of Y back into X to account for
+	;;; subtracting one too many. Knuth shows that the guess is wrong on
+	;;; the order of 3/b, where b is the base (2 to the digit-size power)
+	;;; -- pretty rarely.
+	 (try-bignum-truncate-guess (guess len-y low-x-digit)
+	   (declare (type bignum-index low-x-digit len-y)
+		    (type bignum-element-type guess))
+	   (let ((carry-digit 0)
+		 (borrow 1)
+		 (i low-x-digit))
+	     (declare (type bignum-element-type carry-digit)
+		      (type bignum-index i)
+		      (fixnum borrow))
+	     ;; Multiply guess and divisor, subtracting from dividend
+	     ;; simultaneously.
 	     (dotimes (j len-y)
-	       (multiple-value-bind (v k)
-		   (%add-with-carry (%bignum-ref *truncate-y* j)
-				    (%bignum-ref *truncate-x* i)
-				    carry)
-		 (declare (type bignum-element-type v))
-		 (setf (%bignum-ref *truncate-x* i) v)
-		 (setf carry k))
+	       (multiple-value-bind (high-digit low-digit)
+		   (%multiply-and-add guess
+				      (%bignum-ref truncate-y j)
+				      carry-digit)
+		 (declare (type bignum-element-type high-digit low-digit))
+		 (setf carry-digit high-digit)
+		 (multiple-value-bind (x temp-borrow)
+		     (%subtract-with-borrow (%bignum-ref truncate-x i)
+					    low-digit
+					    borrow)
+		   (declare (type bignum-element-type x)
+			    (fixnum temp-borrow))
+		   (setf (%bignum-ref truncate-x i) x)
+		   (setf borrow temp-borrow)))
 	       (incf i))
-	     (setf (%bignum-ref *truncate-x* i)
-		   (%add-with-carry (%bignum-ref *truncate-x* i) 0 carry)))
-	   (%subtract-with-borrow guess 1 1)))))
+	     (setf (%bignum-ref truncate-x i)
+		   (%subtract-with-borrow (%bignum-ref truncate-x i)
+					  carry-digit borrow))
+	     ;; See whether guess is off by one, adding one
+	     ;; Y back in if necessary.
+	     (cond ((%digit-0-or-plusp (%bignum-ref truncate-x i))
+		    guess)
+		   (t
+		    ;; If subtraction has negative result, add one
+		    ;; divisor value back in. The guess was one too
+		    ;; large in magnitude.
+		    (let ((i low-x-digit)
+			  (carry 0))
+		      (dotimes (j len-y)
+			(multiple-value-bind (v k)
+			    (%add-with-carry (%bignum-ref truncate-y j)
+					     (%bignum-ref truncate-x i)
+					     carry)
+			  (declare (type bignum-element-type v))
+			  (setf (%bignum-ref truncate-x i) v)
+			  (setf carry k))
+			(incf i))
+		      (setf (%bignum-ref truncate-x i)
+			    (%add-with-carry (%bignum-ref truncate-x i)
+					     0 carry)))
+		    (%subtract-with-borrow guess 1 1)))))
+	;;; This returns the amount to shift y to place a one in the
+	;;; second highest bit. Y must be positive. If the last digit
+	;;; of y is zero, then y has a one in the previous digit's
+	;;; sign bit, so we know it will take one less than digit-size
+	;;; to get a one where we want. Otherwise, we count how many
+	;;; right shifts it takes to get zero; subtracting this value
+	;;; from digit-size tells us how many high zeros there are
+	;;; which is one more than the shift amount sought.
+	;;;
+	;;; Note: This is exactly the same as one less than the
+	;;; integer-length of the last digit subtracted from the
+	;;; digit-size.
+	;;;
+	;;; We shift y to make it sufficiently large that doing the
+	;;; 2*digit-size by digit-size %FLOOR calls ensures the quotient and
+	;;; remainder fit in digit-size.
+	 (shift-y-for-truncate (y)
+	   (let* ((len (%bignum-length y))
+		  (last (%bignum-ref y (1- len))))
+	     (declare (type bignum-index len)
+		      (type bignum-element-type last))
+	     (- digit-size (integer-length last) 1)))
+	 ;;; Stores two bignums into the truncation bignum buffers,
+	 ;;; shifting them on the way in. This assumes x and y are
+	 ;;; positive and at least two in length, and it assumes
+	 ;;; truncate-x and truncate-y are one digit longer than x and
+	 ;;; y.
+	 (shift-and-store-truncate-buffers (x len-x y len-y shift)
+	   (declare (type bignum-index len-x len-y)
+		    (type (integer 0 (#.digit-size)) shift))
+	   (cond ((zerop shift)
+		  (bignum-replace truncate-x x :end1 len-x)
+		  (bignum-replace truncate-y y :end1 len-y))
+		 (t
+		  (bignum-ashift-left-unaligned x 0 shift (1+ len-x)
+						truncate-x)
+		  (bignum-ashift-left-unaligned y 0 shift (1+ len-y)
+						truncate-y))))) ;; LABELS
+      ;;; Divide X by Y returning the quotient and remainder. In the
+      ;;; general case, we shift Y to set up for the algorithm, and we
+      ;;; use two buffers to save consing intermediate values. X gets
+      ;;; destructively modified to become the remainder, and we have
+      ;;; to shift it to account for the initial Y shift. After we
+      ;;; multiple bind q and r, we first fix up the signs and then
+      ;;; return the normalized results.
+      (let* ((x-plusp (%bignum-0-or-plusp x (%bignum-length x)))
+	     (y-plusp (%bignum-0-or-plusp y (%bignum-length y)))
+	     (x (if x-plusp x (negate-bignum x nil)))
+	     (y (if y-plusp y (negate-bignum y nil)))
+	     (len-x (%bignum-length x))
+	     (len-y (%bignum-length y)))
+	(multiple-value-bind (q r)
+	    (cond ((< len-y 2)
+		   (bignum-truncate-single-digit x len-x y))
+		  ((plusp (bignum-compare y x))
+		   (let ((res (%allocate-bignum len-x)))
+		     (dotimes (i len-x)
+		       (setf (%bignum-ref res i) (%bignum-ref x i)))
+		     (values 0 res)))
+		  (t
+		   (let ((len-x+1 (1+ len-x)))
+		     (setf truncate-x (%allocate-bignum len-x+1))
+		     (setf truncate-y (%allocate-bignum (1+ len-y)))
+		     (let ((y-shift (shift-y-for-truncate y)))
+		       (shift-and-store-truncate-buffers x len-x y
+							 len-y y-shift)
+		       (values (return-quotient-leaving-remainder len-x+1
+								  len-y)
+			       ;; Now that RETURN-QUOTIENT-LEAVING-REMAINDER
+			       ;; has executed, we just tidy up the remainder
+			       ;; (in TRUNCATE-X) and return it.
+			       (cond
+				 ((zerop y-shift)
+				  (let ((res (%allocate-bignum len-y)))
+				    (declare (type bignum-type res))
+				    (bignum-replace res truncate-x :end2 len-y)
+				    (%normalize-bignum res len-y)))
+				 (t
+				  (shift-right-unaligned
+				   truncate-x 0 y-shift len-y
+				   ((= j res-len-1)
+				    (setf (%bignum-ref res j)
+					  (%ashr (%bignum-ref truncate-x i)
+						 y-shift))
+				    (%normalize-bignum res res-len))
+				   res))))))))
+	  (let ((quotient (cond ((eq x-plusp y-plusp) q)
+				((typep q 'fixnum) (the fixnum (- q)))
+				(t (negate-bignum-in-place q))))
+		(rem (cond (x-plusp r)
+			   ((typep r 'fixnum) (the fixnum (- r)))
+			   (t (negate-bignum-in-place r)))))
+	    (values (if (typep quotient 'fixnum)
+			quotient
+			(%normalize-bignum quotient (%bignum-length quotient)))
+		    (if (typep rem 'fixnum)
+			rem
+			(%normalize-bignum rem (%bignum-length rem))))))))))
 
-;;; This returns a guess for the next division step. Y1 is the highest y
-;;; digit, and y2 is the second to highest y digit. The x... variables are
-;;; the three highest x digits for the next division step.
-;;;
-;;; From Knuth, our guess is either all ones or x-i and x-i-1 divided by y1,
-;;; depending on whether x-i and y1 are the same. We test this guess by
-;;; determining whether guess*y2 is greater than the three high digits of x
-;;; minus guess*y1 shifted left one digit:
-;;;    ------------------------------
-;;;   |    x-i    |   x-i-1  | x-i-2 |
-;;;    ------------------------------
-;;;    ------------------------------
-;;; - | g*y1 high | g*y1 low |   0   |
-;;;    ------------------------------
-;;;		...		  <   guess*y2     ???
-;;; If guess*y2 is greater, then we decrement our guess by one and try again.
-;;; This returns a guess that is either correct or one too large.
-(defun bignum-truncate-guess (y1 y2 x-i x-i-1 x-i-2)
-  (declare (type bignum-element-type y1 y2 x-i x-i-1 x-i-2))
-  (let ((guess (if (%digit-compare x-i y1)
-		   all-ones-digit
-		   (%floor x-i x-i-1 y1))))
-    (declare (type bignum-element-type guess))
-    (loop
-      (multiple-value-bind (high-guess*y1 low-guess*y1) (%multiply guess y1)
-	(declare (type bignum-element-type low-guess*y1 high-guess*y1))
-	(multiple-value-bind (high-guess*y2 low-guess*y2)
-	    (%multiply guess y2)
-	  (declare (type bignum-element-type high-guess*y2 low-guess*y2))
-	  (multiple-value-bind (middle-digit borrow)
-	      (%subtract-with-borrow x-i-1 low-guess*y1 1)
-	    (declare (type bignum-element-type middle-digit)
-		     (fixnum borrow))
-	    ;; Supplying borrow of 1 means there was no borrow, and we know
-	    ;; x-i-2 minus 0 requires no borrow.
-	    (let ((high-digit (%subtract-with-borrow x-i high-guess*y1 borrow)))
-	      (declare (type bignum-element-type high-digit))
-	      (if (and (%digit-compare high-digit 0)
-		       (or (%digit-greater high-guess*y2 middle-digit)
-			   (and (%digit-compare middle-digit high-guess*y2)
-				(%digit-greater low-guess*y2 x-i-2))))
-		  (setf guess (%subtract-with-borrow guess 1 1))
-		  (return guess)))))))))
-
-;;; This returns the amount to shift y to place a one in the second highest
-;;; bit. Y must be positive. If the last digit of y is zero, then y has a
-;;; one in the previous digit's sign bit, so we know it will take one less
-;;; than digit-size to get a one where we want. Otherwise, we count how many
-;;; right shifts it takes to get zero; subtracting this value from digit-size
-;;; tells us how many high zeros there are which is one more than the shift
-;;; amount sought.
-;;;
-;;; Note: This is exactly the same as one less than the integer-length of the
-;;; last digit subtracted from the digit-size.
-;;;
-;;; We shift y to make it sufficiently large that doing the 2*digit-size
-;;; by digit-size %FLOOR calls ensures the quotient and remainder fit in
-;;; digit-size.
-(defun shift-y-for-truncate (y)
-  (let* ((len (%bignum-length y))
-	 (last (%bignum-ref y (1- len))))
-    (declare (type bignum-index len)
-	     (type bignum-element-type last))
-    (- digit-size (integer-length last) 1)))
-
-;;; Stores two bignums into the truncation bignum buffers, shifting them on the
-;;; way in. This assumes x and y are positive and at least two in length, and
-;;; it assumes *truncate-x* and *truncate-y* are one digit longer than x and y.
-(defun shift-and-store-truncate-buffers (x len-x y len-y shift)
-  (declare (type bignum-index len-x len-y)
-	   (type (integer 0 (#.digit-size)) shift))
-  (cond ((zerop shift)
-	 (bignum-replace *truncate-x* x :end1 len-x)
-	 (bignum-replace *truncate-y* y :end1 len-y))
-	(t
-	 (bignum-ashift-left-unaligned x 0 shift (1+ len-x) *truncate-x*)
-	 (bignum-ashift-left-unaligned y 0 shift (1+ len-y) *truncate-y*))))
 
 ;;;; %FLOOR primitive for BIGNUM-TRUNCATE