From d652bd096516a2217e9c273da88ae9c63386a743 Mon Sep 17 00:00:00 2001
From: Raymond Toy <toy.raymond@gmail.com>
Date: Wed, 4 Jul 2018 13:31:21 -0700
Subject: [PATCH] Rename functions to use the new version by default.

Update tests to reflect the change in names.
---
 src/code/bignum.lisp | 63 +++++++++++++++++++++++++++++---------------
 tests/bignum.lisp    | 12 ++++-----
 2 files changed, 48 insertions(+), 27 deletions(-)

diff --git a/src/code/bignum.lisp b/src/code/bignum.lisp
index 1fcaab32a..61eaaf215 100644
--- a/src/code/bignum.lisp
+++ b/src/code/bignum.lisp
@@ -884,7 +884,12 @@ down to individual words.")
 	  (negate-bignum-in-place result))
 	(%normalize-bignum result (1+ (* 2 n))))))))
 
-(defun classical-multiply-bignums (a b)
+;; Bignum multiply using Knuth's algorithm.  We keep this around for
+;; now so we can compare the new algorithm against this to make sure
+;; this are working.
+;;
+;; TODO: Remove this eventually?
+(defun classical-multiply-bignums-knuth (a b)
   (declare (type bignum-type a b))
   (let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
 	 (b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
@@ -916,16 +921,29 @@ down to individual words.")
     (when negate-res (negate-bignum-in-place res))
     (%normalize-bignum res len-res)))
 
-;; Pretend the bignums are actually unsigned, do an unsigned multiply
-;; and then correct the result.  This is based on the algorithm in
-;; Hacker's Delight.
-(defun classical-multiply-bignum-hd (a b)
+;; Classical multiplication of bignums using Knuth's algorithm
+;; modified to handle signed bignums.  Pretend the bignums are
+;; actually unsigned, do an unsigned multiply and then correct the
+;; result.  This is based on the algorithm in Hacker's Delight.
+;;
+;; Let a[n] and b[n] represent the individual bits of each bignum with
+;; M being the number of bits in a and N being the number of bits in
+;; b.  If these are interpreted as an unsigned number, then we are
+;; multiplying numbers
+;;
+;;  (a + 2^M*a[M-1})*(b + 2^N*b[N-1])
+;;    = a*b + 2^M*u[M-1]*b + 2^N*b[N-1]*a + 2^(M+N)*a[M-1]*b[M-1]
+;;
+;; To get the desired result, we need to subtract out the term
+;; 2^M*u[M-1]*b + 2^N*b[N-1]*a from the product.  The last term
+;; doesn't need to subtracted because we know the product fits in M+N
+;; bits and this term is beyond that.
+(defun classical-multiply-bignums (a b)
   (declare (type bignum-type a b))
   (let* ((len-a (%bignum-length a))
 	 (len-b (%bignum-length b))
 	 (len-res (+ len-a len-b))
 	 (res (%allocate-bignum len-res)))
-	 (negate-res (not (eq a-plusp b-plusp))))
     (declare (type bignum-index len-a len-b len-res))
     ;; Unsigned multiply
     (dotimes (i len-a)
@@ -937,36 +955,39 @@ down to individual words.")
 		 (type bignum-element-type carry-digit x))
 	(dotimes (j len-b)
 	  (multiple-value-bind (big-carry res-digit)
-			       (%multiply-and-add x (%bignum-ref b j)
-						  (%bignum-ref res k)
-						  carry-digit)
+	      (%multiply-and-add x (%bignum-ref b j)
+				 (%bignum-ref res k)
+				 carry-digit)
 	    (declare (type bignum-element-type big-carry res-digit))
 	    (setf (%bignum-ref res k) res-digit)
 	    (setf carry-digit big-carry)
 	    (incf k)))
 	(setf (%bignum-ref res k) carry-digit)))
     ;; Apply corrections if either of the arguments is negative.
-    ;; If a < 0, subtract b*2^M from the result
-    (unless(%bignum-0-or-plusp a)
+    (unless (%bignum-0-or-plusp a len-a)
       (let ((borrow 1))
 	(dotimes (j len-b)
+	  (declare (type bignum-index j))
 	  (let ((index (+ j len-a)))
+	    (declare (type bignum-index index))
 	    (multiple-value-bind (d borrow-out)
-		(%subtract-with-borrow (bignum-ref res index)
-				       (bignum-ref b j)
+		(%subtract-with-borrow (%bignum-ref res index)
+				       (%bignum-ref b j)
 				       borrow)
-	      (setf (bignum-ref res index) d)
+	      (setf (%bignum-ref res index) d)
 	      (setf borrow borrow-out))))))
-    (unless (%bignum-0-or-plusp b)
+    (unless (%bignum-0-or-plusp b len-b)
       (let ((borrow 1))
 	(dotimes (j len-a)
+	  (declare (type bignum-index j))
 	  (let ((index (+ j len-b)))
-	  (multiple-value-bind (d borrow-out)
-	      (%subtract-with-borrow (bignum-ref res index)
-				     (bignum-ref a j)
-				     borrow)
-	    (setf (bignum-ref res index) d)
-	    (setf borrow borrow-out)))))
+	    (declare (type bignum-index index))
+	    (multiple-value-bind (d borrow-out)
+		(%subtract-with-borrow (%bignum-ref res index)
+				       (%bignum-ref a j)
+				       borrow)
+	      (setf (%bignum-ref res index) d)
+	      (setf borrow borrow-out))))))
     (%normalize-bignum res len-res)))
 
 (defparameter *min-karatsuba-bits* 512
diff --git a/tests/bignum.lisp b/tests/bignum.lisp
index 783d59d45..46e7ab743 100644
--- a/tests/bignum.lisp
+++ b/tests/bignum.lisp
@@ -21,8 +21,8 @@
       (dotimes (k 100)
 	(let* ((r1 (gen-bignum range (random 2 rng)))
 	       (r2 (gen-bignum range (random 2 rng)))
-	       (prod-knuth (bignum::classical-multiply-bignums r1 r2))
-	       (prod-hd (bignum::classical-multiply-bignum-hd r1 r2)))
+	       (prod-knuth (bignum::classical-multiply-bignums-knuth r1 r2))
+	       (prod-hd (bignum::classical-multiply-bignums r1 r2)))
 	  (assert-equal prod-knuth prod-hd r1 r2))))))
 
 
@@ -45,12 +45,12 @@
 	     (r2 (gen-bignum range 1)) res)
 	(time
 	 (dotimes (k reps)
-	   (declare (fixnum k)) (setf res
-				      (bignum::classical-multiply-bignums r1 r2))))
+	   (declare (fixnum k))
+	   (setf res (bignum::classical-multiply-bignums-knuth r1 r2))))
 	(print res)
 	(time
 	 (dotimes (k reps)
-	   (declare (fixnum k)) (setf res
-				      (bignum::classical-multiply-bignum-hd r1 r2))))
+	   (declare (fixnum k))
+	   (setf res (bignum::classical-multiply-bignums r1 r2))))
 	(print res)))))
 
-- 
GitLab