From b878c8a81a2149a40a865470dbaecb2049c82a7c Mon Sep 17 00:00:00 2001
From: rtoy <rtoy>
Date: Fri, 20 Aug 2004 13:35:03 +0000
Subject: [PATCH] Incorporate Juhu Snellman's modifications to the Karatsuba
 multiplier so that it doesn't have to round up the size of the numbers to the
 next power of 2 size (words).

---
 code/bignum.lisp | 455 +++++++++++++++++++++++++----------------------
 1 file changed, 247 insertions(+), 208 deletions(-)

diff --git a/code/bignum.lisp b/code/bignum.lisp
index d396d1294..13ad14810 100644
--- a/code/bignum.lisp
+++ b/code/bignum.lisp
@@ -5,7 +5,7 @@
 ;;; Carnegie Mellon University, and has been placed in the public domain.
 ;;;
 (ext:file-comment
-  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/bignum.lisp,v 1.37 2004/08/04 20:30:47 rtoy Exp $")
+  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/bignum.lisp,v 1.38 2004/08/20 13:35:03 rtoy Exp $")
 ;;;
 ;;; **********************************************************************
 ;;;
@@ -569,213 +569,249 @@
 classical multiplication algorithm instead of recursing all the way
 down to individual words.")
 
-(defun karatsuba-internal (u v ret d result)
+(defun karatsuba-internal (u v ret n d result)
   (declare (type bignum-index u v ret)
-	   (type bignum-index d)
+	   (type bignum-index d n)
 	   (type bignum-type result)
 	   (optimize (speed 3) (safety 0) (debug 0)))
 
-  (labels ((bignum-multiply-in-place (u v ret d)
-	     ;; Clear out the result area, which might contain garbage for
-	     ;; other recursions!
-	     (do ((i ret (1+ i))
-		  (end (+ ret (* 2 d))))
-		 ((>= i end))
-	       (setf (%bignum-ref result i) 0))
-	     ;; Stolen from bignum.lisp multiply-bignum: a straight forward
-	     ;; multiplication via the classical schoolbook algorithm.
-	     (dotimes (i d)
-	       (declare (type bignum-index i))
-	       (let ((carry-digit 0)
-		     (x (%bignum-ref result (+ u i)))
-		     (k (+ ret i)))
-		 (declare (type bignum-index k)
-			  (type bignum-element-type carry-digit x))
-		 (do ((j 0 (1+ j))
-		      (v+j v (1+ v+j)))
-		     ((>= j d))
-		   (declare (type bignum-index j v+j))
-		   (multiple-value-bind (big-carry res-digit)
-		       (%multiply-and-add x
-					  (%bignum-ref result v+j)
-					  (%bignum-ref result k)
-					  carry-digit)
-		     (declare (type bignum-element-type big-carry res-digit))
-		     (setf (%bignum-ref result k) res-digit)
-		     (setf carry-digit big-carry)
-		     (incf k)))
-		 (setf (%bignum-ref result k) carry-digit))))
-	   (abs-diff (n x y res result)
-	     (declare (type bignum-index n)
-		      (type bignum-index x y res))
-	     (flet ((bigger-p ()
-		      ;; Return T if X >= Y.
-		      (do ((i (1- (+ x n)) (1- i))
-			   (j (1- (+ y n)) (1- j)))
-			  (())
-			(declare (type bignum-index i j))
-			(let ((x-digit (%bignum-ref result i))
-			      (y-digit (%bignum-ref result j)))
-			  (when (/= x-digit y-digit)
-			    (return (>= x-digit y-digit)))))))
-	       (declare (inline bigger-p))
-	       (let ((sign 1))
-		 (when (not (bigger-p))
-		   (rotatef x y)
-		   (setf sign -1))
-		 (let ((borrow 1))
-		   (do ((i x (1+ i))
-			(j y (1+ j))
-			(k res (1+ k))
-			(end (+ x n)))
-		       ((>= i end))
-		     (declare (type bignum-index i j k end))
-		     (multiple-value-bind (diff new-borrow)
-			 (%subtract-with-borrow (%bignum-ref result i)
-						(%bignum-ref result j)
-						borrow)
-		       (setf (%bignum-ref result k) diff)
-		       (setf borrow new-borrow)))
-		   sign))))
-	   (propagate-carry (start end carry)
-	     (declare (type bignum-type result)
-		      (type bignum-index start end)
-		      (type (integer -4 4) carry))
-	     (cond ((minusp carry)
-		    (multiple-value-bind (sum borrow)
-			(%subtract-with-borrow (%bignum-ref result start)
-					       (- carry) 1)
-		      (setf (%bignum-ref result start) sum)
-		      (do ((i (1+ start) (1+ i)))
-			  ((or (>= i end) (= borrow 1)) (- borrow 1))
-			(declare (type bignum-index i))
-			(multiple-value-bind (sum next-borrow)
-			    (%subtract-with-borrow (%bignum-ref result i)
-						   0 borrow)
-			  (setf borrow next-borrow)
-			  (setf (%bignum-ref result i) sum)))))
-		   (t
-		    (multiple-value-bind (sum carry)
-			(%add-with-carry (%bignum-ref result start) carry 0)
-		      (setf (%bignum-ref result start) sum)
-		      (do ((i (1+ start) (1+ i)))
-			  ((or (>= i end) (zerop carry)) carry)
-			(declare (type bignum-index i))
-			(multiple-value-bind (sum next-carry)
-			    (%add-with-carry (%bignum-ref result i) 0 carry)
-			  (setf carry next-carry)
-			  (setf (%bignum-ref result i) sum)))))))
-	   (stageC-add ()
-	     (let ((carry11 0)
-		   (carry12 0)
-		   (carry2 0)
-		   (d/2 (ash d -1)))
-	       (do ((i0 ret (1+ i0))
-		    (i1 (+ ret d/2) (1+ i1))
-		    (i2 (+ ret d) (1+ i2))
-		    (i3 (+ ret d d/2) (1+ i3))
-		    (i4 (+ ret d d) (1+ i4))
-		    (end (+ ret d/2)))
-		   ((>= i0 end))
-		 (declare (type bignum-index i0 i1 i2 i3 i4))
-		 (let ((b (%bignum-ref result i1)))
-		   (multiple-value-bind (sum c1)
-		       (%add-with-carry b (%bignum-ref result i0) carry11)
-		     (multiple-value-bind (sum c2)
-			 (%add-with-carry sum (%bignum-ref result i3) carry12)
-		       (setf carry11 c1)
-		       (setf carry12 c2)
-		       (setf (%bignum-ref result i1) sum)))
-		   (multiple-value-bind (sum carry)
-		       (%add-with-carry b (%bignum-ref result i4) carry2)
-		     (setf (%bignum-ref result i2) sum)
-		     (setf carry2 carry))))
-	       (+ (propagate-carry (+ ret d) (+ ret d d/2) (+ carry11 carry12))
-		  carry2)))
-	   (stageC-subtract ()
-	     (let ((carry11 0)
-		   (borrow12 1)
-		   (borrow2 1)
-		   (d/2 (ash d -1)))
-	       (do ((i0 ret (1+ i0))
-		    (i1 (+ ret d/2) (1+ i1))
-		    (i2 (+ ret d) (1+ i2))
-		    (i3 (+ ret d d/2) (1+ i3))
-		    (i4 (+ ret d d) (1+ i4))
-		    (end (+ ret d/2)))
-		   ((>= i0 end))
-		 (declare (type bignum-index i0 i1 i2 i3 i4))
-		 (let ((b (%bignum-ref result i1)))
-		   (multiple-value-bind (sum c1)
-		       (%add-with-carry b (%bignum-ref result i0) carry11)
-		     (multiple-value-bind (sum b2)
-			 (%subtract-with-borrow sum (%bignum-ref result i3)
-						borrow12)
-		       (setf carry11 c1)
-		       (setf borrow12 b2)
-		       (setf (%bignum-ref result i1) sum)))
-		   (multiple-value-bind (sum borrow)
-		       (%subtract-with-borrow b (%bignum-ref result i4)
-					      borrow2)
-		     (setf (%bignum-ref result i2) sum)
-		     (setf borrow2 borrow))))
-	       (+ (propagate-carry (+ ret d) (+ ret d d/2)
-				   (+ carry11 borrow12 -1))
-		  (- borrow2 1))))
-	   (stageE-add (stageC-carry)
-	     ;; Sum the above pieces, result to R0124.
-	     (let ((carry1 0)
-		   (carry21 0)
-		   (carry22 0)
-		   (d/2 (ash d -1)))
-	       (declare (type (mod 2) carry1 carry21 carry22))
-	       (do ((i0 ret (1+ i0))
-		    (i1 (+ ret d/2) (1+ i1))
-		    (i2 (+ ret d) (1+ i2))
-		    (i3 (+ ret d d/2) (1+ i3))
-		    (i4 (+ ret d d) (1+ i4))
-		    (end (+ ret d/2)))
-		   ((>= i0 end))
-		 (declare (type bignum-index i0 i1 i2 i3 i4))
-		 (multiple-value-bind (sum c1)
-		     (%add-with-carry (%bignum-ref result i1)
-				      (%bignum-ref result i3)
-				      carry1)
-		   (setf carry1 c1)
-		   (setf (%bignum-ref result i1) sum))
-		 (let ((e (%bignum-ref result i4)))
+  (let* ((odd (logand d 1))
+	 (origd d)
+	 (d/2- (ash d -1))
+	 (d/2 (+ d/2- odd))
+	 (d (* d/2 2))
+	 (n/2 (ash n -1)))
+    (labels ((bignum-multiply-in-place (u v ret d result cleard)
+	       (declare (type bignum-index u v ret)
+			(type bignum-index d)
+			(type bignum-type result)
+			(optimize (speed 3) (safety 0) (debug 0)))
+	       ;; Clear out the result area, which might contain garbage for
+	       ;; other recursions!
+	       (do ((i ret (1+ i))
+		    (end (+ ret (* 2 cleard))))
+		   ((>= i end))
+		 (setf (%bignum-ref result i) 0))
+	       ;; Stolen from bignum.lisp multiply-bignum: a straight forward
+	       ;; multiplication via the classical schoolbook algorithm.
+	       (dotimes (i d)
+		 (declare (type bignum-index i))
+		 (let ((carry-digit 0)
+		       (x (%bignum-ref result (+ u i)))
+		       (k (+ ret i)))
+		   (declare (type bignum-index k)
+			    (type bignum-element-type carry-digit x))
+		   (do ((j d (1- j))
+			(v+j v (1+ v+j)))
+		       ((zerop j))
+		     (declare (type bignum-index j v+j))
+		     (multiple-value-bind (big-carry res-digit)
+			 (%multiply-and-add x
+					    (%bignum-ref result v+j)
+					    (%bignum-ref result k)
+					    carry-digit)
+		       (declare (type bignum-element-type big-carry res-digit))
+		       (setf (%bignum-ref result k) res-digit)
+		       (setf carry-digit big-carry)
+		       (incf k)))
+		   (setf (%bignum-ref result k) carry-digit))))
+	     (abs-diff (n1 x n2 y res result)
+	       (declare (type bignum-index n1 n2)
+			(type bignum-index x y res))
+	       (flet ((bigger-p ()
+			;; We guarantee in the caller that n1 is either
+			;; n2 or (1+ n2)
+			(when (> n1 n2)
+			  ;;(assert (= n1 (1+ n2)))
+			  (if (zerop (%bignum-ref result (+ x n1 -1)))
+			      (progn
+				(setf (%bignum-ref result (+ res n1 -1)) 0)
+				(decf n1))				
+			      (return-from bigger-p t)))
+			;; Return T if X >= Y.		      
+			(do ((i (1- (+ x n1)) (1- i))
+			     (j (1- (+ y n2)) (1- j)))
+			    (())
+			  (declare (type bignum-index i j))
+			  (let ((x-digit (%bignum-ref result i))
+				(y-digit (%bignum-ref result j)))
+			    (when (/= x-digit y-digit)
+			      (return-from bigger-p (>= x-digit y-digit))))
+			  ;; Need to test for the exit condition here
+			  ;; instead of in the DO headers because when
+			  ;; X or Y is 0 and the numbers are equal I
+			  ;; or J will be decremented to -1, which
+			  ;; isn't of type bignum-index.
+			  (when (= i x)
+			    (return-from bigger-p t)))))
+		 (declare (inline bigger-p))
+		 (let ((sign 1))
+		   (when (not (bigger-p))		     
+		     (rotatef x y)
+		     (rotatef n1 n2)
+		     (setf sign -1))
+		   (let ((borrow 1))
+		     (do ((i x (1+ i))
+			  (j y (1+ j))
+			  (k res (1+ k))
+			  (end (+ x n1))
+			  (end2 (+ y n2)))
+			 ((>= i end))
+		       (declare (type bignum-index i j k end))
+		       (multiple-value-bind (diff new-borrow)
+			   (%subtract-with-borrow (%bignum-ref result i)
+						  (if (>= j end2)
+						      0
+						      (%bignum-ref result j))
+						  borrow)
+			 (setf (%bignum-ref result k) diff)
+			 (setf borrow new-borrow)))
+		     sign))))
+	     (propagate-carry (start end carry)
+	       (declare (type bignum-type result)
+			(type bignum-index start end)
+			(type (integer -4 4) carry))
+	       (cond ((minusp carry)
+		      (multiple-value-bind (sum borrow)
+			  (%subtract-with-borrow (%bignum-ref result start)
+						 (- carry) 1)
+			(setf (%bignum-ref result start) sum)
+			(do ((i (1+ start) (1+ i)))
+			    ((or (>= i end) (= borrow 1)) (- borrow 1))
+			  (declare (type bignum-index i))
+			  (multiple-value-bind (sum next-borrow)
+			      (%subtract-with-borrow (%bignum-ref result i)
+						     0 borrow)
+			    (setf borrow next-borrow)
+			    (setf (%bignum-ref result i) sum)))))
+		     (t
+		      (multiple-value-bind (sum carry)
+			  (%add-with-carry (%bignum-ref result start) carry 0)
+			(setf (%bignum-ref result start) sum)
+			(do ((i (1+ start) (1+ i)))
+			    ((or (>= i end) (zerop carry)) carry)
+			  (declare (type bignum-index i))
+			  (multiple-value-bind (sum next-carry)
+			      (%add-with-carry (%bignum-ref result i) 0 carry)
+			    (setf carry next-carry)
+			    (setf (%bignum-ref result i) sum)))))))
+	     (stageC-add ()
+	       (let ((carry11 0)
+		     (carry12 0)
+		     (carry2 0))
+		 (do ((i0 ret (1+ i0))
+		      (i1 (+ ret d/2) (1+ i1))
+		      (i2 (+ ret d) (1+ i2))
+		      (i3 (+ ret n/2) (1+ i3))
+		      (i4 (+ ret n/2 d/2) (1+ i4))
+		      (end (+ ret d/2)))
+		     ((>= i0 end))
+		   (declare (type bignum-index i0 i1 i2 i3 i4))
+		   (let ((b (%bignum-ref result i1)))
+		     (multiple-value-bind (sum c1)
+			 (%add-with-carry b (%bignum-ref result i0) carry11)
+		       (multiple-value-bind (sum c2)
+			   (%add-with-carry sum (%bignum-ref result i3) carry12)
+			 (setf carry11 c1)
+			 (setf carry12 c2)
+			 (setf (%bignum-ref result i1) sum)))
+		     (multiple-value-bind (sum carry)
+			 (%add-with-carry b (%bignum-ref result i4) carry2)
+		       (setf (%bignum-ref result i2) sum)
+		       (setf carry2 carry))))
+		 (+ (propagate-carry (+ ret d) (+ ret n/2) (+ carry11 carry12))
+		    carry2)))
+	     (stageC-subtract ()
+	       (let ((carry11 0)
+		     (borrow12 1)
+		     (borrow2 1))
+		 (do ((i0 ret (1+ i0))
+		      (i1 (+ ret d/2) (1+ i1))
+		      (i2 (+ ret d) (1+ i2))
+		      (i3 (+ ret n/2) (1+ i3))
+		      (i4 (+ ret n/2 d/2) (1+ i4))
+		      (end (+ ret d/2)))
+		     ((>= i0 end))
+		   (declare (type bignum-index i0 i1 i2 i3 i4))
+		   (let ((b (%bignum-ref result i1)))
+		     (multiple-value-bind (sum c1)
+			 (%add-with-carry b (%bignum-ref result i0) carry11)
+		       (multiple-value-bind (sum b2)
+			   (%subtract-with-borrow sum (%bignum-ref result i3)
+						  borrow12)
+			 (setf carry11 c1)
+			 (setf borrow12 b2)
+			 (setf (%bignum-ref result i1) sum)))
+		     (multiple-value-bind (sum borrow)
+			 (%subtract-with-borrow b (%bignum-ref result i4)
+						borrow2)
+		       (setf (%bignum-ref result i2) sum)
+		       (setf borrow2 borrow))))
+		 (+ (propagate-carry (+ ret d) (+ ret n/2)
+				     (+ carry11 borrow12 -1))
+		    (- borrow2 1))))
+	     (stageE-add (stageC-carry)
+	       ;; Sum the above pieces, result to R0124.
+	       (let ((carry1 0)
+		     (carry21 0)
+		     (carry22 0))
+		 (declare (type (mod 2) carry1 carry21 carry22))
+		 (do ((i0 ret (1+ i0))
+		      (i1 (+ ret d/2) (1+ i1))
+		      (i2 (+ ret d) (1+ i2))
+		      (i3 (+ ret d d/2) (1+ i3))
+		      (i4 (+ ret d d) (1+ i4))
+		      (end (+ ret d/2)))
+		     ((>= i0 end))
+		   (declare (type bignum-index i0 i1 i2 i3 i4))
 		   (multiple-value-bind (sum c1)
-		       (%add-with-carry e (%bignum-ref result i2) carry21)
-		     (multiple-value-bind (sum c2)
-			 (%add-with-carry sum (%bignum-ref result i3) carry22)
-		       (setf carry21 c1)
-		       (setf carry22 c2)
-		       (setf (%bignum-ref result i2) sum)))
-		   (setf (%bignum-ref result i3) e)))
-	       (assert (zerop (propagate-carry (+ ret d) (+ ret d d) carry1)))
-	       (assert (zerop (propagate-carry (+ ret d d/2) (+ ret d d)
-					       (+ carry21 carry22 stageC-carry)))))))
-
-    (cond ((<= d *karatsuba-classical-cutoff*)
-	   (bignum-multiply-in-place u v ret d))
-	  (t
-	   (let* ((d/2 (ash d -1))
-		  (ret1 (+ ret d/2))
-		  (ret3 (+ ret d d/2))
-		  (u-hi (+ u d/2))
-		  (v-hi (+ v d/2))
-		  (diff-sign (* (abs-diff d/2 u u-hi ret result)
-				(abs-diff d/2 v v-hi ret1 result))))
-	     (declare (type bignum-index d/2 ret1 ret3 u-hi v-hi)
-		      (type (member -1 1) diff-sign))
-	     (karatsuba-internal ret ret1 ret3 d/2 result)
-	     (karatsuba-internal u v ret d/2 result)
-	     (let ((stagec-carry (if (minusp diff-sign)
-				     (stagec-add)
-				     (stagec-subtract))))
-	       (karatsuba-internal u-hi v-hi ret3 d/2 result)
-	       (stageE-add stageC-carry))))))
-  result)
+		       (%add-with-carry (%bignum-ref result i1)
+					(%bignum-ref result i3)
+					carry1)
+		     (setf carry1 c1)
+		     (setf (%bignum-ref result i1) sum))
+		   (let ((e (%bignum-ref result i4)))
+		     (multiple-value-bind (sum c1)
+			 (%add-with-carry e (%bignum-ref result i2) carry21)
+		       (multiple-value-bind (sum c2)
+			   (%add-with-carry sum (%bignum-ref result i3) carry22)
+			 (setf carry21 c1)
+			 (setf carry22 c2)
+			 (setf (%bignum-ref result i2) sum)))
+		     (setf (%bignum-ref result i3) e)))
+		 (propagate-carry (+ ret d) (+ ret d d) carry1)
+		 (propagate-carry (+ ret d d/2)
+						 (+ ret d d)
+						 (+ carry21 carry22 stageC-carry)))))
+      (cond ((<= d (ext:truly-the fixnum *karatsuba-classical-cutoff*))
+	     (bignum-multiply-in-place u v ret origd result d))
+	    (t
+	     (let* ((ret1 (+ ret d/2))
+		    (ret3 (+ ret n/2))
+		    (u-hi (+ u d/2))
+		    (v-hi (+ v d/2))
+		    (diff-sign (* (abs-diff d/2 u d/2- u-hi ret result)
+				  (abs-diff d/2 v d/2- v-hi ret1 result))))
+	       (declare (type bignum-index d/2 ret1 ret3 u-hi v-hi)
+			(type (member -1 1) diff-sign))
+	       (do ((i ret3 (1+ i)))
+		   ((>= i (+ ret3 d)))
+		 (setf (%bignum-ref result i) 0))
+	       (karatsuba-internal ret ret1 ret3 n/2 d/2 result)
+	       (do ((i ret (1+ i)))
+		   ((>= i ret3))
+		 (setf (%bignum-ref result i) 0))	       
+	       (karatsuba-internal u v ret n/2 d/2 result)
+	       (let ((stagec-carry (if (minusp diff-sign)
+				       (stagec-add)
+				       (stagec-subtract))))
+		 (setf ret3 (+ ret d d/2))
+		 (karatsuba-internal u-hi v-hi ret3 n/2 d/2- result)
+		 (when (/= d/2 d/2-)
+		   (setf (%bignum-ref result (+ ret3 (* d/2- 2))) 0)
+		   (setf (%bignum-ref result (+ ret3 (* d/2- 2) 1)) 0))
+		 (stageE-add stageC-carry))))))
+    result))
 
 
 ;;; karatsuba  --  Interface.
@@ -789,7 +825,7 @@ down to individual words.")
 ;;;
 (defun karatsuba (x y)
   (declare (type bignum-type x y)
-	   (optimize (speed 3) (safety 0) (debug 2)))
+	   (optimize (speed 3) (safety 0) (debug 3)))
   (flet ((power-of-two (n)
 	   ;; Compute the smallest power of two greater than or equal
 	   ;; to the given number.
@@ -806,7 +842,8 @@ down to individual words.")
 	   (negate-res (not (eq x-plusp y-plusp)))
 	   (x-words (ceiling (integer-length x) 32))
 	   (y-words (ceiling (integer-length y) 32))
-	   (n (power-of-two (max x-words y-words))))
+	   (d (max x-words y-words))
+	   (n (power-of-two d)))
       (declare (type bignum-index n))
 
       (let ((result (%allocate-bignum (* 5 n))))
@@ -825,13 +862,15 @@ down to individual words.")
 	(setf (%bignum-ref result save)
 	      (%bignum-ref y k)))
       ;; Do it.
-      (karatsuba-internal (* 3 n) (* 4 n) 0 n result)
+      (karatsuba-internal (* 3 n) (* 4 n) 0 (* n 3)
+			  (max x-words y-words)
+			  result)
       ;; Convert the result into a bignum.
       (let ((new-len (1+ (* 2 n)))
 	    (n5 (* n 5)))
-	(do ((k (* 2 n) (1+ k)))
+	(do ((k (* 2 d) (1+ k)))
 	    ((>= k n5))
-	  (assert (< k (%bignum-length result)))
+	  ;;(assert (< k (%bignum-length result)))
 	  (setf (%bignum-ref result k) 0))
 	(%bignum-set-length result new-len)
 	(when negate-res
-- 
GitLab