Skip to content
Snippets Groups Projects
Commit b878c8a8 authored by rtoy's avatar rtoy
Browse files

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).
parent dd281a02
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment