| ... |
... |
@@ -884,7 +884,12 @@ down to individual words.") |
|
884
|
884
|
(negate-bignum-in-place result))
|
|
885
|
885
|
(%normalize-bignum result (1+ (* 2 n))))))))
|
|
886
|
886
|
|
|
887
|
|
-(defun classical-multiply-bignums (a b)
|
|
|
887
|
+;; Bignum multiply using Knuth's algorithm. We keep this around for
|
|
|
888
|
+;; now so we can compare the new algorithm against this to make sure
|
|
|
889
|
+;; this are working.
|
|
|
890
|
+;;
|
|
|
891
|
+;; TODO: Remove this eventually?
|
|
|
892
|
+(defun classical-multiply-bignums-knuth (a b)
|
|
888
|
893
|
(declare (type bignum-type a b))
|
|
889
|
894
|
(let* ((a-plusp (%bignum-0-or-plusp a (%bignum-length a)))
|
|
890
|
895
|
(b-plusp (%bignum-0-or-plusp b (%bignum-length b)))
|
| ... |
... |
@@ -916,16 +921,29 @@ down to individual words.") |
|
916
|
921
|
(when negate-res (negate-bignum-in-place res))
|
|
917
|
922
|
(%normalize-bignum res len-res)))
|
|
918
|
923
|
|
|
919
|
|
-;; Pretend the bignums are actually unsigned, do an unsigned multiply
|
|
920
|
|
-;; and then correct the result. This is based on the algorithm in
|
|
921
|
|
-;; Hacker's Delight.
|
|
922
|
|
-(defun classical-multiply-bignum-hd (a b)
|
|
|
924
|
+;; Classical multiplication of bignums using Knuth's algorithm
|
|
|
925
|
+;; modified to handle signed bignums. Pretend the bignums are
|
|
|
926
|
+;; actually unsigned, do an unsigned multiply and then correct the
|
|
|
927
|
+;; result. This is based on the algorithm in Hacker's Delight.
|
|
|
928
|
+;;
|
|
|
929
|
+;; Let a[n] and b[n] represent the individual bits of each bignum with
|
|
|
930
|
+;; M being the number of bits in a and N being the number of bits in
|
|
|
931
|
+;; b. If these are interpreted as an unsigned number, then we are
|
|
|
932
|
+;; multiplying numbers
|
|
|
933
|
+;;
|
|
|
934
|
+;; (a + 2^M*a[M-1})*(b + 2^N*b[N-1])
|
|
|
935
|
+;; = a*b + 2^M*u[M-1]*b + 2^N*b[N-1]*a + 2^(M+N)*a[M-1]*b[M-1]
|
|
|
936
|
+;;
|
|
|
937
|
+;; To get the desired result, we need to subtract out the term
|
|
|
938
|
+;; 2^M*u[M-1]*b + 2^N*b[N-1]*a from the product. The last term
|
|
|
939
|
+;; doesn't need to subtracted because we know the product fits in M+N
|
|
|
940
|
+;; bits and this term is beyond that.
|
|
|
941
|
+(defun classical-multiply-bignums (a b)
|
|
923
|
942
|
(declare (type bignum-type a b))
|
|
924
|
943
|
(let* ((len-a (%bignum-length a))
|
|
925
|
944
|
(len-b (%bignum-length b))
|
|
926
|
945
|
(len-res (+ len-a len-b))
|
|
927
|
946
|
(res (%allocate-bignum len-res)))
|
|
928
|
|
- (negate-res (not (eq a-plusp b-plusp))))
|
|
929
|
947
|
(declare (type bignum-index len-a len-b len-res))
|
|
930
|
948
|
;; Unsigned multiply
|
|
931
|
949
|
(dotimes (i len-a)
|
| ... |
... |
@@ -937,36 +955,39 @@ down to individual words.") |
|
937
|
955
|
(type bignum-element-type carry-digit x))
|
|
938
|
956
|
(dotimes (j len-b)
|
|
939
|
957
|
(multiple-value-bind (big-carry res-digit)
|
|
940
|
|
- (%multiply-and-add x (%bignum-ref b j)
|
|
941
|
|
- (%bignum-ref res k)
|
|
942
|
|
- carry-digit)
|
|
|
958
|
+ (%multiply-and-add x (%bignum-ref b j)
|
|
|
959
|
+ (%bignum-ref res k)
|
|
|
960
|
+ carry-digit)
|
|
943
|
961
|
(declare (type bignum-element-type big-carry res-digit))
|
|
944
|
962
|
(setf (%bignum-ref res k) res-digit)
|
|
945
|
963
|
(setf carry-digit big-carry)
|
|
946
|
964
|
(incf k)))
|
|
947
|
965
|
(setf (%bignum-ref res k) carry-digit)))
|
|
948
|
966
|
;; Apply corrections if either of the arguments is negative.
|
|
949
|
|
- ;; If a < 0, subtract b*2^M from the result
|
|
950
|
|
- (unless(%bignum-0-or-plusp a)
|
|
|
967
|
+ (unless (%bignum-0-or-plusp a len-a)
|
|
951
|
968
|
(let ((borrow 1))
|
|
952
|
969
|
(dotimes (j len-b)
|
|
|
970
|
+ (declare (type bignum-index j))
|
|
953
|
971
|
(let ((index (+ j len-a)))
|
|
|
972
|
+ (declare (type bignum-index index))
|
|
954
|
973
|
(multiple-value-bind (d borrow-out)
|
|
955
|
|
- (%subtract-with-borrow (bignum-ref res index)
|
|
956
|
|
- (bignum-ref b j)
|
|
|
974
|
+ (%subtract-with-borrow (%bignum-ref res index)
|
|
|
975
|
+ (%bignum-ref b j)
|
|
957
|
976
|
borrow)
|
|
958
|
|
- (setf (bignum-ref res index) d)
|
|
|
977
|
+ (setf (%bignum-ref res index) d)
|
|
959
|
978
|
(setf borrow borrow-out))))))
|
|
960
|
|
- (unless (%bignum-0-or-plusp b)
|
|
|
979
|
+ (unless (%bignum-0-or-plusp b len-b)
|
|
961
|
980
|
(let ((borrow 1))
|
|
962
|
981
|
(dotimes (j len-a)
|
|
|
982
|
+ (declare (type bignum-index j))
|
|
963
|
983
|
(let ((index (+ j len-b)))
|
|
964
|
|
- (multiple-value-bind (d borrow-out)
|
|
965
|
|
- (%subtract-with-borrow (bignum-ref res index)
|
|
966
|
|
- (bignum-ref a j)
|
|
967
|
|
- borrow)
|
|
968
|
|
- (setf (bignum-ref res index) d)
|
|
969
|
|
- (setf borrow borrow-out)))))
|
|
|
984
|
+ (declare (type bignum-index index))
|
|
|
985
|
+ (multiple-value-bind (d borrow-out)
|
|
|
986
|
+ (%subtract-with-borrow (%bignum-ref res index)
|
|
|
987
|
+ (%bignum-ref a j)
|
|
|
988
|
+ borrow)
|
|
|
989
|
+ (setf (%bignum-ref res index) d)
|
|
|
990
|
+ (setf borrow borrow-out))))))
|
|
970
|
991
|
(%normalize-bignum res len-res)))
|
|
971
|
992
|
|
|
972
|
993
|
(defparameter *min-karatsuba-bits* 512
|