GitLab

at 2018-07-04T09:54:08-07:00

Bignum multiply without consing temp space

The current bignum multiplier creates temp space to hold the absolute
value of the bignums and then negates the result (in-place) at the
end.

Instead, use the algorithm from Hacker's Delight that pretends the
numbers are unsigned, does the unsigned multiply and finally corrects
the result.  No extra memory is needed for this.

@@ -3,6 +3,7 @@ variables:
    version: "2018-03-x86"
  linux-runner:
 +  image: ubuntu:14.04
    tags:
      - linux
    before_script:

@@ -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,6 +921,73 @@ down to individual words.")
      (when negate-res (negate-bignum-in-place res))
      (%normalize-bignum res len-res)))
 +;; 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)))
 +    (declare (type bignum-index len-a len-b len-res))
 +    ;; Unsigned multiply
 +    (dotimes (i len-a)
 +      (declare (type bignum-index i))
 +      (let ((carry-digit 0)
 +	    (x (%bignum-ref a i))
 +	    (k i))
 +	(declare (type bignum-index k)
 +		 (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)
 +	    (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)))
 +    (flet ((apply-correction (neg-arg neg-len pos-arg pos-len)
 +	     ;; Applies the correction by basically subtracting out
 +	     ;; 2^M*b where M is the length (in bits) of b and b is
 +	     ;; the positive term in pos-arg.  neg-arg is the negative
 +	     ;; arg.
 +	     (let ((borrow 1))
 +	       (dotimes (j pos-len)
 +		 (declare (type bignum-index j))
 +		 (let ((index (+ j neg-len)))
 +		   (declare (type bignum-index index))
 +		   (multiple-value-bind (d borrow-out)
 +		       (%subtract-with-borrow (%bignum-ref res index)
 +					      (%bignum-ref pos-arg j)
 +					      borrow)
 +		     (setf (%bignum-ref res index) d)
 +		     (setf borrow borrow-out)))))))
 +      ;; Apply corrections if either of the arguments is negative.
 +      (unless (%bignum-0-or-plusp a len-a)
 +	;; A is negative
 +	(apply-correction a len-a b len-b))
 +      (unless (%bignum-0-or-plusp b len-b)
 +	;; B is negative
 +	(apply-correction b len-b a len-a)))
 +    (%normalize-bignum res len-res)))
++
  (defparameter *min-karatsuba-bits* 512
    "Use Karatsuba if the bignums have at least this many bits")

 +;;; Tests for the bignum operations
++
 +(defpackage :bignum-tests
 +  (:use :cl :lisp-unit))
++
 +(in-package #:bignum-tests)
++
 +(define-test hd-mult.same-size
 +  "Test bignum multiplier"
 +  (:tag :bignum-tests)
 +  ;; x and y are randomly generated 128 integers. No particular reason
 +  ;; for these values, except that they're bignums.
 +  (let ((x 248090201001762284446997112921270181259)
 +	(y 313102667534462314033767199170708979663)
 +	(prod 77677703722812705876871716049945873590003455155145426220435549433670954735717))
 +    ;; Verify the we get the right results for various signed values of x and y.
 +    (assert-equal prod (* x y))
 +    (assert-equal (- prod) (* (- x) y))
 +    (assert-equal (- prod) (* x (- y)))
 +    (assert-equal prod (* (- x) (- y)))
 +    ;; Nake sure it's commutative
 +    (assert-equal prod (* y x))
 +    (assert-equal (- prod) (* y (- x)))
 +    (assert-equal (- prod) (* (- y) x))
 +    (assert-equal prod (* (- y) (- x)))))
++
 +(define-test hd-mult.diff-size
 +  "Test bignum multiplier"
 +  (:tag :bignum-tests)
 +  ;; x is a randomly generated bignum.  y is a small bignum.
 +  (let ((x 248090201001762284446997112921270181259)
 +	(y (1+ most-positive-fixnum))
 +	(prod 133192412470079431258262755675409306410924638208))
 +    ;; Verify the we get the right results for various signed values of x and y.
 +    (assert-equal prod (* x y))
 +    (assert-equal (- prod) (* (- x) y))
 +    (assert-equal (- prod) (* x (- y)))
 +    (assert-equal prod (* (- x) (- y)))
 +    ;; Nake sure it's commutative
 +    (assert-equal prod (* y x))
 +    (assert-equal (- prod) (* y (- x)))
 +    (assert-equal (- prod) (* (- y) x))
 +    (assert-equal prod (* (- y) (- x)))))
++
++
 +(define-test hd-mult.random
 +  "Test bignum multiplier with random values"
 +  (:tag :bignum-tests)
 +  (let ((rng (kernel::make-random-object :state (kernel:init-random-state)
 +					 :rand 0
 +					 :cached-p nil))
 +	(range (ash 1 128)))
 +    (flet ((gen-bignum (x sign)
 +	     (do ((r (random x rng) (random x rng)))
 +		 ((typep r 'bignum)
 +		  (if (zerop sign)
 +		      r
 +		      (- r))))))
 +      (dotimes (k 100)
 +	(let* ((r1 (gen-bignum range (random 2 rng)))
 +	       (r2 (gen-bignum range (random 2 rng)))
 +	       (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))))))
++
++
 +;; Just for simple timing tests so we can redo the timing tests if needed.
 +#+nil
 +(define-test hd-timing
 +  "Test execution time"
 +  (:tag :bignum-tests)
 +  (let ((rng (kernel::make-random-object :state
 +					 (kernel:init-random-state)
 +					 :rand 0 :cached-p nil))
 +	(range (ash 1 128))
 +	(reps 10000))
 +    (flet ((gen-bignum (x sign)
 +	     (do ((r (random x rng) (random x rng)))
 +		 ((typep r 'bignum)
 +		  (if (zerop sign)
 +		      r (- r))))))
 +      (let* ((r1 (gen-bignum range 1))
 +	     (r2 (gen-bignum range 1)) res)
 +	(time
 +	 (dotimes (k reps)
 +	   (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-bignums r1 r2))))
 +	(print res)))))
++

@@ -397,6 +397,10 @@
        (sleep 5)
        (assert-eql :exited (ext:process-status p)))))
 +;; For some reason this used to work linux CI but not doesn't.  But
 +;; this test passes on my Fedora and debian systesm.
 +;; See issue #64.
 +#-linux
  (define-test issue.41.1
      (:tag :issues)
    (issue-41-tester unix:sigstop))

Raymond Toy pushed to branch master at cmucl / cmucl

Commits:

4 changed files:

Changes: