Better bignum multiplication

Do not cons positive bignums at the start of the loop, the result can
be adjusted in place by subtraction at the end.
On arm64 and x86-64, write the whole loop in assembly, not just the
first iteration.

src/code/bignum.lisp

@@ -485,53 +485,65 @@
                         %normalize-bignum-buffer))
 
 ;;;; multiplication
+(declaim (inline %subtract-bignum-in-place))
+(defun %subtract-bignum-in-place (a b offset len)
+  (declare (type bignum a b)
+           (type bignum-length offset len)
+           (optimize speed (safety 0)))
+  (let ((borrow 1))
+    (declare (type (integer 0 1) borrow))
+    (dotimes (i len)
+      (declare (type bignum-index i))
+      (let ((a-idx (the bignum-index (+ offset i))))
+        (multiple-value-bind (diff new-borrow)
+            (%subtract-with-borrow (%bignum-ref a a-idx)
+                                   (%bignum-ref b i)
+                                   borrow)
+          (setf (%bignum-ref a a-idx) diff)
+          (setf borrow new-borrow))))))
 
+;;; See x86-64-vm.lisp for variants with mulx
 (defun multiply-bignums (a b)
   (declare (type bignum a b)
            (optimize speed (safety 0)))
-  (let* ((a-plusp (bignum-plus-p a))
+  (let* ((len-a (%bignum-length a))
-         (b-plusp (bignum-plus-p b))
-         (a (if a-plusp a (negate-bignum-not-fully-normalized a)))
-         (b (if b-plusp b (negate-bignum-not-fully-normalized b)))
-         (len-a (%bignum-length a))
          (len-b (%bignum-length b))
          (len-res (+ len-a len-b))
-         (res (%allocate-bignum len-res))
+         (res (%allocate-bignum len-res)))
-         (negate-res (not (eq a-plusp b-plusp))))
-    (declare (type bignum-length len-a len-b len-res))
     (when (> len-a len-b)
       (rotatef a b)
       (rotatef len-a len-b))
+    (sb-c::if-vop-existsp (:named sb-vm::bignum-multiply-loop)
+      (sb-sys:%primitive sb-vm::bignum-multiply-loop a len-a b len-b res)
 
-    ;; The partial result is zero on the first iteration,
+      ;; The partial result is zero on the first iteration,
-    ;; so don't include it. And no need to zero when allocating it.
+      ;; so don't include it. And no need to zero when allocating it.
-    (let ((x (%bignum-ref a 0)))
+      (let ((x (%bignum-ref a 0)))
-      (sb-c::if-vop-existsp (:named sb-vm::bignum-mult-and-add-word-loop)
-        (sb-sys:%primitive sb-vm::bignum-mult-and-add-word-loop b x len-b res)
         (let ((carry-digit 0))
           (declare (fixnum carry-digit))
           (dotimes (index len-b)
             (declare (type bignum-index index))
-            (setf (values carry-digit
+            (setf (values carry-digit (%bignum-ref res index))
-                          (%bignum-ref res index))
                   (%multiply-and-add (%bignum-ref b index) x carry-digit)))
-          (setf (%bignum-ref res len-b) carry-digit))))
+          (setf (%bignum-ref res len-b) carry-digit))
-
+        (loop for i of-type bignum-index from 1 below len-a
-    (loop for i of-type bignum-index from 1 below len-a
+              do
-          do
+              (let ((x (%bignum-ref a i)))
-          (let ((x (%bignum-ref a i))
+                (let ((k i)
-                (k i)
+                      (carry-digit 0))
-                (carry-digit 0))
+                  (declare (type bignum-index k))
-            (declare (type bignum-index k))
+                  (dotimes (j len-b)
-            (dotimes (j len-b)
+                    (setf (values carry-digit (%bignum-ref res k))
-              (setf (values carry-digit (%bignum-ref res k))
+                          (%multiply-and-add x
-                    (%multiply-and-add x
+                                             (%bignum-ref b j)
-                                       (%bignum-ref b j)
+                                             (%bignum-ref res k)
-                                       (%bignum-ref res k)
+                                             carry-digit))
-                                       carry-digit))
+                    (incf k))
-              (incf k))
+                  (setf (%bignum-ref res k) carry-digit))))))
-            (setf (%bignum-ref res k) carry-digit)))
+    (unless (bignum-plus-p a)
-    (when negate-res (negate-bignum-in-place res))
+      (%subtract-bignum-in-place res b len-a len-b))
+    (unless (bignum-plus-p b)
+      (%subtract-bignum-in-place res a len-b len-a))
     (%normalize-bignum res len-res)))
 
 (defun multiply-bignum-and-fixnum (bignum fixnum)
@@ -542,25 +554,26 @@
         ((eql fixnum -1)
          (- bignum))
         (t
-         (let* ((bignum-plus-p (bignum-plus-p bignum))
+         (let* ((bignum-len (%bignum-length bignum))
-                (fixnum-plus-p (not (minusp fixnum)))
+                (abs-fixnum (abs fixnum))
-                (bignum (if bignum-plus-p bignum (negate-bignum-not-fully-normalized bignum)))
-                (bignum-len (%bignum-length bignum))
-                (fixnum (if fixnum-plus-p fixnum (- fixnum)))
                 (result (%allocate-bignum (1+ bignum-len))))
            (declare (type bignum bignum result)
-                    (type bignum-element-type fixnum))
+                    (type bignum-element-type abs-fixnum))
            (sb-c::if-vop-existsp (:named sb-vm::bignum-mult-and-add-word-loop)
-             (sb-sys:%primitive sb-vm::bignum-mult-and-add-word-loop bignum fixnum bignum-len result)
+             (sb-sys:%primitive sb-vm::bignum-mult-and-add-word-loop bignum abs-fixnum bignum-len result)
              (let ((carry-digit 0))
                (declare (fixnum carry-digit))
                (dotimes (index bignum-len)
                  (declare (type bignum-index index))
                  (setf (values carry-digit
                                (%bignum-ref result index))
-                       (%multiply-and-add (%bignum-ref bignum index) fixnum carry-digit)))
+                       (%multiply-and-add (%bignum-ref bignum index) abs-fixnum carry-digit)))
                (setf (%bignum-ref result bignum-len) carry-digit)))
-           (unless (eq bignum-plus-p fixnum-plus-p)
+           (unless (bignum-plus-p bignum)
+             (setf (%bignum-ref result bignum-len)
+                   (logand (- (%bignum-ref result bignum-len) abs-fixnum)
+                           most-positive-word)))
+           (when (minusp fixnum)
              (negate-bignum-in-place result))
            (%normalize-bignum result (1+ bignum-len))))))
 

src/code/x86-64-vm.lisp

@@ -352,6 +352,7 @@
 
 (in-package :sb-bignum)
 
+#+x86-64
 (sb-vm::def-variant multiply-bignum-and-fixnum :bmi2 (bignum fixnum)
   (declare (type bignum bignum) (type fixnum fixnum)
            (optimize speed (safety 0)))
@@ -360,17 +361,35 @@
         ((eql fixnum -1)
          (- bignum))
         (t
-         (let* ((bignum-plus-p (bignum-plus-p bignum))
+         (let* ((bignum-len (%bignum-length bignum))
-                (fixnum-plus-p (not (minusp fixnum)))
+                (abs-fixnum (abs fixnum))
-                (bignum (if bignum-plus-p bignum (negate-bignum-not-fully-normalized bignum)))
-                (bignum-len (%bignum-length bignum))
-                (fixnum (if fixnum-plus-p fixnum (- fixnum)))
                 (result (%allocate-bignum (1+ bignum-len))))
            (declare (type bignum bignum result)
-                    (type bignum-element-type fixnum))
+                    (type bignum-element-type abs-fixnum))
-           (sb-sys:%primitive sb-vm::bignum-mulx-and-add-word-loop bignum fixnum bignum-len result)
+           (sb-sys:%primitive sb-vm::bignum-mulx-and-add-word-loop bignum abs-fixnum bignum-len result)
-           (unless (eq bignum-plus-p fixnum-plus-p)
+           (unless (bignum-plus-p bignum)
+             (setf (%bignum-ref result bignum-len)
+                   (logand (- (%bignum-ref result bignum-len) abs-fixnum)
+                           most-positive-word)))
+           (when (minusp fixnum)
              (negate-bignum-in-place result))
            (%normalize-bignum result (1+ bignum-len))))))
 
+(sb-vm::def-variant multiply-bignums :bmi2 (a b)
+  (declare (type bignum a b)
+           (optimize speed (safety 0)))
+  (let* ((len-a (%bignum-length a))
+         (len-b (%bignum-length b))
+         (res (%allocate-bignum (+ len-a len-b))))
+    (when (> len-a len-b)
+      (rotatef a b)
+      (rotatef len-a len-b))
+
+    (sb-sys:%primitive sb-vm::bignum-mulx-multiply-loop a len-a b len-b res)
+    (unless (bignum-plus-p a)
+      (%subtract-bignum-in-place res b len-a len-b))
+    (unless (bignum-plus-p b)
+      (%subtract-bignum-in-place res a len-b len-a))
+    (%normalize-bignum res (+ len-a len-b))))
+
 (in-package :sb-vm)

src/compiler/arm64/arith.lisp

@@ -2036,6 +2036,77 @@
     (inst adc hi hi zr-tn)
     (inst str hi (@ r (extend index :lsl 3)))))
 
+(define-vop (bignum-multiply-loop)
+  (:args (a* :scs (descriptor-reg))
+         (la :scs (unsigned-reg))
+         (b* :scs (descriptor-reg))
+         (lb :scs (unsigned-reg))
+         (r* :scs (descriptor-reg)))
+  (:arg-types bignum unsigned-num bignum unsigned-num bignum)
+  (:temporary (:sc descriptor-reg) a b-base b r-row-base r)
+  (:temporary (:sc unsigned-reg) outer-len inner-len
+                                 x b-digit r-digit lo hi carry)
+  (:generator 40
+    (inst add-sub a a* (- (* bignum-digits-offset n-word-bytes) other-pointer-lowtag))
+    (inst add-sub b-base b* (- (* bignum-digits-offset n-word-bytes) other-pointer-lowtag))
+    (inst add-sub r-row-base r* (- (* bignum-digits-offset n-word-bytes) other-pointer-lowtag))
+
+
+    (inst ldr x (@ a n-word-bytes :post-index))
+    (inst mov b b-base)
+    (inst mov r r-row-base)
+    (inst mov inner-len lb)
+    (inst adds hi zr-tn zr-tn) ;; clear carry
+
+    ROW0-LOOP
+    (inst ldr b-digit (@ b n-word-bytes :post-index))
+    (inst mul lo b-digit x)
+    (inst adcs lo lo hi)
+    (inst umulh hi b-digit x)
+    (inst str lo (@ r n-word-bytes :post-index))
+    (inst sub inner-len inner-len 1)
+    (inst cbnz inner-len ROW0-LOOP)
+
+    (inst adc hi hi zr-tn)
+    (inst str hi (@ r))
+
+    (inst subs outer-len la 1)
+    (inst b :eq DONE)
+
+    (inst add r-row-base r-row-base n-word-bytes)
+
+    OUTER-LOOP
+    (inst ldr x (@ a n-word-bytes :post-index))
+    (inst mov b b-base)
+    (inst mov r r-row-base)
+    (inst mov inner-len lb)
+    (inst adds carry zr-tn zr-tn) ;; clear carry
+
+    INNER-LOOP
+    (inst ldr b-digit (@ b n-word-bytes :post-index))
+    (inst ldr r-digit (@ r))
+
+    (inst mul lo b-digit x)
+    (inst umulh hi b-digit x)
+
+    (inst adcs lo lo carry)
+    (inst adc carry hi zr-tn)
+
+    (inst adds lo lo r-digit)
+
+    (inst str lo (@ r n-word-bytes :post-index))
+    (inst sub inner-len inner-len 1)
+    (inst cbnz inner-len INNER-LOOP)
+
+    (inst adc carry carry zr-tn)
+    (inst str carry (@ r))
+
+    (inst add r-row-base r-row-base n-word-bytes)
+    (inst subs outer-len outer-len 1)
+    (inst b :ne OUTER-LOOP)
+
+    DONE))
+
 (define-vop (bignum-mult-and-add-3-arg)
   (:translate sb-bignum:%multiply-and-add)
   (:policy :fast-safe)

src/compiler/x86-64/arith.lisp

@@ -4158,6 +4158,74 @@
     (inst adc hi 0)
     (inst mov (ea #1# r index 8) hi)))
 
+(define-vop (bignum-mulx-multiply-loop)
+  (:args (a* :scs (descriptor-reg))
+         (la :scs (unsigned-reg unsigned-stack))
+         (b* :scs (descriptor-reg))
+         (lb :scs (unsigned-reg unsigned-stack))
+         (r* :scs (descriptor-reg)))
+  (:arg-types bignum unsigned-num bignum unsigned-num bignum)
+  (:temporary (:sc unsigned-reg :offset rdx-offset) rdx)
+  (:temporary (:sc unsigned-reg) a b-end r-end index carry lo hi outer-len)
+  (:generator 40
+    (inst mov :dword index lb)
+    (inst lea b-end (ea #1=(- (* bignum-digits-offset n-word-bytes) other-pointer-lowtag) b* index 8))
+    (inst lea r-end (ea #1# r* index 8))
+
+    (inst mov rdx (ea #1# a*))
+
+    (inst neg index)
+    (zeroize carry) ;; clears cf
+
+    ROW0-LOOP
+    (inst mulx hi lo (ea 0 b-end index 8))
+    (inst adc lo carry)
+    (inst mov carry hi)
+    (inst mov (ea 0 r-end index 8) lo)
+
+    (inst inc index)
+    (inst jmp :nz ROW0-LOOP)
+
+    (inst adc carry 0)
+    (inst mov (ea r-end) carry)
+
+    (inst mov :dword index la)
+    (inst dec :dword index)
+    (inst jmp :z DONE)
+    (inst mov :dword outer-len index)
+
+    (inst lea a (ea (+ #1# n-word-bytes) a*))
+    (inst add r-end n-word-bytes)
+
+    OUTER-LOOP
+    (inst mov rdx (ea a))
+    (inst add a n-word-bytes)
+
+    (inst mov :dword index lb)
+    (inst neg index)
+    (zeroize carry) ;; clears cf
+
+    INNER-LOOP
+    (inst mulx hi lo (ea 0 b-end index 8))
+
+    (inst adc lo carry)
+    (inst adc hi 0)
+
+    (inst add (ea 0 r-end index 8) lo)
+    (inst mov carry hi)
+
+    (inst inc index)
+    (inst jmp :nz INNER-LOOP)
+
+    (inst adc carry 0)
+    (inst mov (ea r-end) carry)
+
+    (inst add r-end n-word-bytes)
+    (inst dec :dword outer-len)
+    (inst jmp :nz OUTER-LOOP)
+
+    DONE))
+
 (define-vop (bignum-mult-and-add-3-arg)
   (:translate sb-bignum:%multiply-and-add)
   (:policy :fast-safe)