Raymond Toy pushed to branch rtoy-bignum-mult-less-consing at cmucl / cmucl
Commits: d652bd09 by Raymond Toy at 2018-07-04T13:31:21-07:00 Rename functions to use the new version by default.
Update tests to reflect the change in names.
- - - - -
2 changed files:
- src/code/bignum.lisp - tests/bignum.lisp
Changes:
===================================== src/code/bignum.lisp ===================================== --- 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
===================================== tests/bignum.lisp ===================================== --- 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)))))
View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/commit/d652bd096516a2217e9c273da8...